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ABSTRACT 


Sensor nodes in a wireless sensor network (WSN) can establish a link with a 
UAV by using beamforming techniques to form a random array with position errors. The 
position errors’ effect in the array performance is examined using a MATLAB-based 


simulation model. 


In order to spread the processing and communication load among the nodes, two 
new distributed algorithms for beamforming in WSN, based on the least squares (LS) 
approximation of the desired array response, are proposed. The first is a distributed 
implementation of the QR decomposition, and the second is an iterative method for 
solving the LS problem. Results indicate that the processing load is effectively shared 
among the nodes. Especially, in the second approach, the processing load can be lower 
than that of the centralized approach, depending on the algorithm’s convergence. For 
both algorithms, the tradeoff for the ability to spread the processing load is the increased 
communication cost, which could cause an overall increase in the total power 
consumption in the network. However, the average power per participating sensor node is 
still lower than that required by the cluster head in the centralized approach. 
Consequently, the network’s susceptibility to failures due to excessive power 


consumption is greatly reduced. 
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EXECUTIVE SUMMARY 


A wireless sensor network (WSN) consists of a large number of microsensors, 
each having limited battery lifetime and restricted communication and computing 
capabilities. Recent advances in the integrated circuit technology have allowed the 
production of lightweight and inexpensive sensor nodes, which have a range of 
capabilities, such as sensing, processing and communication. WSNs can have many 
applications in both commercial and military environments. A set of sensor nodes, which 
can be deployed easily and quickly by an unmanned aerial vehicle (UAV), for example, 
can be used for monitoring the battlefield environment, sensing for a wide range of 
targets, especially in the case where the area of interest is inaccessible or there is high 


risk of human loss. 


Once deployed, the sensor nodes can collect the desired information and transmit 
it to the UAV. Although a single sensor node cannot transmit its data directly to the VAV 
due to the limited range of coverage, several of them can coordinate their transmissions 
in order to form an array and thus substantially increase the range of coverage. Since the 
sensors are randomly deployed, it is unlikely that they form topologies that permit the 
formation of arrays with equally spaced elements. As a result, there are position errors 


with respect to an array with equally spaced elements. 


A simulation model was developed in the MATLAB environment to analyze the 
effect of these position errors on the array performance. The simulations showed that the 
sidelobe levels in the array response increase as a function of the error in the element 
location. Specifically, as the mean deviation from the ideal position was increased, the 
average sidelobe magnitudes also increased. These results are in agreement with the 
theoretical analysis of random arrays, found in the literature. The degradation of the array 
performance can be largely eliminated using a Least Squares (LS) beamformer, which 
computes the weights that best approximate a given desired response. This beamformer 
can also efficiently suppress potential interfering signals coming from directions other 


than the signal’s. 


XVil 


Since reliability and robustness in a sensor network environment are desired, the 
processing load must be effectively spread among the sensor nodes. Centralized 
approaches assign the entire processing load to a single node whereas in a distributed 
approach, the processing tasks are split into smaller processes, which are then allocated to 
the participating sensor nodes. Two fully distributed approaches to beamforming in WSN 
were presented in this work, and they are both based on the LS approximation of the 
desired response. The first is a distributed implementation of the QR decomposition, and 


the second is an iterative method of computing the weights in the LS sense. 


The performance of the distributed methods was compared to the centralized LS 
approach using the processing and communication costs as metrics. The results indicate 
that the processing load is effectively shared among the nodes. Especially in the second 
method, the processing load is a function of the algorithm’s convergence and can be 
lower compared to that of the centralized approach, subject to the speed of convergence. 
For both algorithms, the tradeoff for the ability to spread the processing load is the 
increased communication cost, which could cause an overall increase in the total power 
consumption in the network. This total power, however, is shared among the sensor 
nodes; therefore, the average power expended by a participating sensor node in the 
distributed implementation is lower than the power required by the cluster head in the 
centralized approach. Consequently, the network’s susceptibility to failures due to 


excessive power consumption is greatly reduced. 


XVill 


I. INTRODUCTION 


A. INTRODUCTION TO WIRELESS SENSOR NETWORKS 


A wireless sensor network (WSN) consists of a large number of microsensors, 
each having limited battery lifetime and, therefore, restricted communication and 
computing capabilities [1], [2] . Recent advances in the integrated circuit technology have 
allowed the production of lightweight and inexpensive sensor nodes, which have a range 
of capabilities, such as sensing, processing and communication. If they are properly 
networked and programmed, these sensor nodes can cooperate in order to perform 


complex signal processing functions [3], [4]. 


The main issue for a WSN is to prolong its operational lifetime as much as 
possible, taking into account the sensors’ power consumption requirements. Stringent 
energy limitations are also a crucial factor when designing signal processing algorithms 
for a WSN [5]. Since such energy restrictions are not taken into account by the signal 
processing methods that are already used in applications other than WSNs, existing 
techniques should be modified in order to conform to the sensor nodes’ specific 
characteristics. Therefore, a major challenge in recent research is the design of signal 
processing and networking operations, which optimize the tradeoff between energy 


efficiency, simplicity, and overall performance, [3]. 


Because microsensors are becoming cheaper and more capable, WSNs will find 
more applications in both commercial and military environments [1], [2]. Future tactical 
operations will involve the deployment of large-scale WSNs in which hundreds or 
thousands of disposable sensor nodes will cooperate in order to achieve the mission 
objective [3]. These nodes can be deployed easily and quickly by an unmanned aerial 
vehicle (UAV), for example, as in Figure 1, which minimizes the risk of human loss. 
Then they can be used for monitoring the battlefield environment, sensing for a wide 
range of targets, such as biological, radioactive, nuclear, chemical and other materials [1]. 
Once deployed, the sensor nodes can collect the desired information and disseminate it to 


a relay node, such as a UAV. Furthermore, taking into account that the sensing 


environment may be harsh and inaccessible for deploying wired networks, there is 


obvious need for developing WSNs consisting of small and disposable sensors. 
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Figure 1. WSN deployed over an area of interest and UAV collecting the desired 
information. 


B. RELATED WORK IN BEAMFORMING AND WIRELESS SENSOR 
NETWORKS 


After forming an ad hoc network and collecting the required data about the target 
of interest, the sensor nodes must establish communication with a UAV, so the acquired 
information can be transmitted to the UAV. Although promising, today’s technology still 
imposes strict limits on the processing and communication capabilities of the sensor 
nodes [6], [7]. Single sensor nodes do not have sufficient power to communicate with an 
overflying UAV. Since the UAV may be required to fly at a high altitude due to the 
hostile nature of the operative environment, the objective of transmitting the collected 


information to the UAV becomes more difficult. 


Although a single sensor node cannot transmit its data directly to the UAV due to 
the limited batter power, several of them can cooperate in order to function as a large 


antenna array and thus substantially increase the transmission range and the data rate. 
2 


This process of combining the signal from different antenna elements in order to form a 
single output of the sensor array is known as beamforming. It has been proven by many 
studies that when an antenna array is properly configured, it can improve the channel 
capacity and extend the range of coverage [8], [9]. It can also reduce the multipath fading 
and the bit error rate (BER); therefore, it results in more reliable communication [8]. 
Additionally, beamforming can adaptively steer the antenna beam towards the UAV, thus 
aiming the radiated energy in the desired direction [10]. Furthermore, the antenna gain is 
proportional to the number of the antenna elements, so the main beam peak power 
density can be of several orders of magnitude higher than that of a single sensor [9]. 
Another useful characteristic of an antenna array is that it can be used in order to perform 
spatio-temporal filtering, thus suppressing potential interference signals coming from 
directions other than the desired direction [11], [12]. In summary, taking the above 
mentioned advantages into account, beamforming in WSNs can meet the objective of 


establishing an efficient communication link between a WSN and a UAV. 


Several algorithms for beamforming exist in the literature [13], [11], [14] and 
many of them are successfully implemented in conventional antenna arrays. [10]. 
Nevertheless, these algorithms for the computation of the weights for the array elements 
cannot be directly implemented in WSNs since there are significant differences between 
WSNs and conventional arrays. For instance, the phased arrays used in RADAR are 
installed permanently on site [15]; thus, the positions of the antenna elements are fixed. 
On the other hand, the sensor nodes in a WSN are usually randomly deployed and their 
relative positions are not predetermined. Due to this random deployment, there are 
position errors, which cause performance deterioration of the antenna array, compared to 
that of an array of equally spaced elements. Moreover, the sensors are prone to frequent 
failures due to limited battery life or due to their vulnerability to environmental 
conditions. Therefore, the topology of the sensor array can change substantially as new 


nodes are added or withdrawn. 


Another significant problem is that, in conventional arrays, where power is not a 
major issue, the beamforming operation is performed in a single processor [16]. All 
necessary information is collected in a central processing node, which is responsible for 
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solving the beamforming problem. However, in a sensor network environment, this is 
neither reliable nor desirable since a single node would be assigned this computationally 
demanding task. Additionally, such centralized implementations create a single point of 
failure, which in turn creates a serious system vulnerability. If this node fails, then the 
beamforming problem has to be solved from the beginning. Thus, the processing load or 
consequently the power usage should be effectively distributed across the sensor network 
[17]. However, an optimum set of participating sensors in the array has to be defined 
since the communication cost for organizing the sensor nodes into an array is prohibitive 


after a certain critical number of nodes [18]. 


Cc. THESIS OBJECTIVE 


The objective of this research is to implement several distributed beamforming 
algorithms, and to evaluate their performance. Throughout this work, the operational 
scenario of Figure | is adopted where several sensor nodes try to communicate with a 
UAV. The beamforming process is not performed in a central node (cluster head), but it 
is split into smaller processes, which then can be allocated to the sensor nodes. The main 
concept is that the processing and communication cost must be shared among the nodes, 
so there is no single point of failure and that the energy of the nodes 1s efficiently used, 
thus extending the WSN lifetime. This work is focused on investigating beamforming 
schemes that have the same performance as well established centralized approaches yet 
offer the advantage of implementation in a distributed fashion thus, increasing the 


network’s robustness and overall performance. 


D. PROPOSED APPROACH TO DISTRIBUTED BEAMFORMING 


Starting from a centralized approach to the Least Squares (LS) solution of the 
beamforming problem where the beamformer is designed in such way that the desired 
array performance is best approximated, two fully distributed methods are proposed. The 


first one is a distributed implementation of the QR decomposition with Householder 
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transformations and derives the exact solution for the array weight vector. The second 
scheme is an iterative method for solving the LS problem, implemented in a distributed 


fashion. 


In order to examine the proposed techniques, a performance analysis of the 
communication and computational costs is developed. These two costs are closely 
connected to the power consumption and provide a reliable test for the algorithms’ 
effectiveness. The resulting array response is examined and compared with the desired 
response. The results from these two distributed implementations are encouraging and 
indicate that they can provide a realistic solution for the beamforming problem in sensor 


networks. 


E. THESIS OUTLINE 


Chapter II introduces the fundamental concepts of the antenna arrays, including a 
description of the uniform linear and planar array. This is followed by an analysis of the 
effects of position errors on the performance of the antenna array and a simulation, which 
confirms the theoretical results. Beamforming in wireless sensor networks is also 
presented along with a specific operational communication scenario which uses a UAV. 
Finally, the framework for evaluating the algorithms’ performance based on the factors 
that affect power consumption, such as the processing and communication cost, is 


developed. 


Chapter III presents two centralized beamforming approaches and evaluates their 
performance. Their advantages and disadvantages are discussed, and their ability to 
mitigate position errors is analyzed using a simulation model developed in a MATLAB 


environment. 


In Chapter IV, two distributed algorithms for beamforming in wireless sensor 
networks are proposed. Their performance is analyzed in terms of efficiently sharing the 
processing and communication cost among the nodes, and the simulation results are 


compared to those obtained by the centralized approaches. 


Chapter V summarizes the significant results of this thesis and provides some 


ideas for extending this work in the future. 


Finally, the Appendix includes the MATLAB code used in the simulations. 


HW. BEAMFORMING IN SENSOR ARRAYS 


In this chapter, the main concepts of antenna arrays are discussed; specifically, the 
uniform linear and planar arrays are presented as well as their array response 
(beampattern). The uniform array is followed by an analysis of the random array with 
position errors. The effects of the antenna element position errors in the main lobe and 
the sidelobe power gain are presented under various assumptions about the statistical 
characteristics of the error. Beamformer implementations are described, particularly the 
narrowband beamformer. Lastly, there is a discussion about communication and 
computational cost in sensor networks as a function of power consumption, which 


indicates the need for distributed algorithms in beamforming. 


A. UNIFORM LINEAR AND PLANAR ARRAYS 


The fundamental concepts of the uniform linear and planar arrays are presented in 
this section along with the basic formulation of beamforming which will be used 


throughout this work. 


tL, One-Dimensional Array 


In Figure 2, a linear array is depicted with M identical, equally spaced elements. 
The spacing between consecutive array elements or sensor nodes in the case of a sensor 
atray is assumed equal to a half wavelength, i.e, d=4/2, where the elements are 
isotropic, meaning their beampattern is omni-directional. Furthermore, it is assumed that 
the array is located far enough from the signal source (e.g., a UAV); thus, it is considered 
to lie in the source’s far field. The array axis is assumed to be in the x-axis. The plane 


wave s(t) arrives at the array at an angle @, with respect to the x-axis (array axis), and 


for this reason, the m+1" element senses the signal earlier than the m™ element. If x,, is 


‘h 


the distance between the m” node and the reference node, which is located at the origin 


of the coordinate system, then the signal s(t) arrives at the m™ element earlier by t,, 


seconds with respect to the reference element. The time difference, which depends on the 
arrival angle 6, or Angle of Arrival (AOA), and the element’s distance from the 


reference point is given by [19]: 


Ln (4,, ) = (1) 
where c is the speed of light. 


Z-axis 


source 
e s(t) 


ee : 





x-axis 


Xm 





Figure 2. An M x1 linear array of equally spaced isotropic elements. 


Each array element is weighted by a complex weight w,, for m=0,1,...,M, 


which multiplies the incoming signal. Adding all the elements’ weighted inputs gives the 


spatial response of the array or array factor F(@,) for any arbitrary angle 6): 


M 
RO)=yiWern (2) 
m=1 


where w. =1.e7" and I. and e””"* are the magnitude and the phase of the 
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complex weights, respectively. Using the wavenumber { = 27/4 and the expression for 


JBXn COS Gay 


the time difference ¢,,, the weights can be written as w,, =I,e and the array 


m? 


factor as 
M FS 
HG) = Dale ee (3) 
m=1 


The weights w,, are carefully selected in order to give the maximum value of the 
array response F(@,) at the desired direction @, and to suppress potential interference 
signals arriving from other directions. Indeed at 6, =6,,, the array response reaches its 


maximum value 


nu i Bx, S j BX, COS ‘a 
FO) be eg ey eM. at Ia vin (4) 
m=l 


m=1 
The set of weights w,, form the weight vector 
w=[w Ww, .. w, | 
while the steering vector or direction vector is defined as 


d(@,) = [l ei Px cos6, eh cos6, Be eh cos 6, id 


and incorporates the location information of the array. Therefore, the array response can 


be expressed as 
F(O,)=w" d(6,). (5) 


Note that the main concept of beamforming is the use of the weights w in order to point 
the array beam towards any desired direction. So, if the desired transmission direction is 


@,,, then the beamformer should set its weights to be 


[PX COS O,,, 
Wee Lert ee (6) 


The beampattern of a 10x1 uniform linear array is shown in Figure 3 where the 


normalized power gain G is defined as [9] 


2 





G(o,) = |F(O,) 


aE or (7) 
max |F(,) 


and it is plotted as a function of the direction @, (in degrees). The array elements are 


identical and isotropic, and the spacing has a fixed value of 2/2 and beam pointing 


angle 6, =0°. The maximum sidelobe is equal to -13 dB and the Half Power Beamwidth 


(HPBW) is about 10.2°. 
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Figure 3. Normalized Power Gain (dB) — Beampattern of a M =10 element array with 


isotropic elements, fixed spacing 2/2 and 6, =0°. 


In general, the sidelobe level decreases with increasing number of elements M 
and approaches the value of -13.3 dB. The HPBW in any plane containing the array axis. 
is given [9] by 


A 
8, 43 = 0.866 — 
3dB Md 
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where Mdl A, 1e., it is valid for long arrays. Therefore the HPBW is a function of 


M and decreases as the number of array elements increase. 


2. Two-dimensional (planar) Array 

The previous discussion of a linear array can be easily expanded to a planar array, 
and similar expressions can be derived for the array response and the power gain. In 
Figure 4, a uniformly spaced planar array is depicted with M x N identical and isotropic 
elements. The spacing between the array elements is assumed equal to half wavelength, 
d=A/2, in both directions while it is assumed again that the array is located in the 


source’s far field. The plane wave s(t) arrives at the array at polar angle @, with respect 
to the z-axis and an azimuth angle ¢ with respect to the x-axis; thus, the (m,n) element 
receives the signal earlier by f,,, seconds compared to the reference element at the origin. 
This time difference ¢,, depends on the angles 6,, ¢ and the element’s position 


(Xnn> Ynn) 1 the array, and is given by [19] 


Xn SIN B, COS P + Yn, SIN A sin gy 


tn (Gy >A) = (8) 


Cc 


Adding all the elements’ weighted inputs gives the array response of the planar 


array F(0,@) for any arbitrary choice of angles @ and @ : 


N M 
= * — {B( Xm SiN OcosP+y,,, Sin O sin p) 
FO, ~) = > y Win’ (9) 
n=| m=1 
where 
_ JB (Xn Si. COS Py + Yin Si Fp Sin Py ) 
Wan _ Lon (10) 


are the complex weights for each element (m,n). In matrix form, the above expression 


can be written as 
F(0,¢)=w' d(0,¢) (11) 


where w isan MN x1 weight vector and 


11 


1 
eh sin Ocos¢+ y,7 sin@sin¢) 
eh sin@ cos + y,3 sin@sing) 
d(0,) = (12) 


e/hGim sin Ocos¢+ y,y sinOsing) 





e/P Cun sinO cos $+ yyy sin@sing) 





is an MN x] steering vector. 


The weights w,,, are again selected in order to give the maximum value of the 
array response F(0,@) in the desired direction(@,¢,). At 0=@, and g=4¢,, the array 


response reaches its maximum value 


F(@,Q) =MN if 7, =1,Vm,n (13) 


where in (13) it is assumed that the excitation of elements (weighting) is uniform. 
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Figure 4. Signal wavefront arriving from direction (@,,@,) at a M x N planar array. 
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In Figure 5, the three-dimensional beampattern of a 10x10 uniform planar array 
is shown, where the normalized power gain G in dB is plotted as a function of the polar 


angle @ and the azimuth angle ¢. A cross section of this 3-D beampattern is plotted in 
Figure 6 where the azimuth angle is constant at g, =45° (direction of arrival), and the 


beampattern varies with the polar angle @. 


Normalized Power Gain (dB) 





¢ (degrees) -100  -100 


0 (degrees) 


Figure 5. Normalized Power Gain (dB): 3-D Beampattern of a 10x10 planar array with 
isotropic elements, equally spaced by 4/2 and direction of signal 
arrival (@, = 30°, d, =45°). 
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Figure 6. Normalized Power Gain (dB): Beampattern of a 10x10 planar array with 
isotropic elements, equally spaced by 4/2 and direction of signal 
arrival (9, = 0°, @ =45°). Azimuth angle is fixed at ¢ =45°. 


The maximum sidelobe is equal to -26 dB, and the HPBW is about 10.2°. In 
general, the sidelobe levels decrease with an increasing number of elements MN and 
approaches the value of -26.6 dB. Similarly, the HPBW decreases continuously as the 


number of elements increase. 


B. RANDOM ARRAY: POSITION ERRORS 


In a random deployment of a sensor array where the sensor nodes are dropped 
randomly over an area of interest, it would be unrealistic to expect formations of 
perfectly spaced planar arrays. An illustrating example comes from the fact that for an 


operating frequency f,. = 1 GHz, the ideal distance between the sensor nodes is 


A/2=15 cm, and even a small displacement of 3 cm yields a position error of 20%. 


Therefore, the performance of a sensor array should be studied using the theoretical 
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analysis of random arrays. In general, the beampattern of a randomly deployed sensor 
array will be affected by position errors, amplitude and phase errors, quantization errors 
and failures of the nodes, which cause a change in the array topology. Throughout this 
work, the main emphasis will be given to the effect of position errors and solutions to this 


problem. 


1. Position and Phase Errors 


There are several references [19], [20], [21], [22] in the literature that deal with 
position errors in random arrays. The random antenna elements’ misplacement causes 
phase errors and mismatches, which yield degraded performance for the array. In [19], an 
analysis of the radiation pattern of a random array with both amplitude and phase errors 
is presented. For an M x N two-dimensional array, assuming that there are no amplitude 


errors, 1.e., the weights w,,, have the same magnitude J = /,,,, the expected increase A, in 


mn? 


the sidelobe level with respect to the main lobe is given by [19] 


ae 
a 2(¢ ~—1) 14 
val oe 


where the phase errors follow a Gaussian distribution with zero mean and variance o, ae 


For example, ina 5x5 planar array, A, is equal to 0.00518 for o,,, =20° and equal to 


0.025 for o,, =40°. Therefore, doubling of the phase error will cause an increase of 


almost 6 dB in the sidelobe level with respect to the main lobe. However, it can be seen 
from (14) that the effect of phase errors can be mitigated by increasing the number of 


array elements in the array. Indeed, A, can be decreased by a factor of two if the number 


of antenna elements is doubled. 
The fractional loss in the main lobe gain due to phase errors is given by [19] 


Siege (15) 
Gy 


where G and G_, are the main lobe power gains with and without the presence of phase 


errors, respectively. Since for a 3 dB reduction in the main lobe level, a phase error 
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standard deviation of o,,=47° is needed, it can be concluded that the main lobe is not 


significantly affected by random position errors. 


2. Random Array Implementation 


In general, the effects of misplacement errors depend strongly on the assumptions 
about the random characteristics of the position deviations from the uniform array. In the 
previous section, the results were derived based on the assumption of Gaussian phase 
errors. However, there are several references in the literature which consider different 
deployments of sensor arrays and consequently different assumptions about the random 
distribution of the phase errors. One such example is included in [22] where the location 
of each node in the sensor array is chosen randomly by a uniform distribution within a 


disk. 


Throughout this work, the position errors will be modeled as a deviation from the 
ideal position of a uniform array. The position errors will be modeled as a uniformly 
distributed random variable with a minimum value of 0 and a maximum value of 
axA/2 where a is the maximum percentage error; therefore, the mean error will be 
axA/4. It is also assumed that there is displacement in both x and y directions as 
indicated in Figure 7. The array axis for randomly placed elements is the best line fit of 
the sensor nodes’ positions and it is assumed to be also the x-axis of the coordinates 
system. In Figure 7 the best line fit of the nodes location is found and it is assumed to be 
the array axis. Then the deviations are defined using this array axis as a reference. The 


beampattern is computed in a plane which is perpendicular to the x-y plane and at ¢, (in 


this case @, = 45°) with respect to the x-axis. 


In Figure 8, the effect of position errors in a 10x1 linear antenna array topology 
is depicted. The sidelobe levels have been increased and consequently the array 
performance has been deteriorated. The mean beampattern of a 10x1 linear array with 


mean position deviation equal to 0.32/2 in both x and y directions is compared to the 
beampattern of the ideal uniform 10x1 array for polar angle 6, =30° and azimuth angle 


g@ =45°. The array axis is in the x-direction and the beampattern is in a plane 
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perpendicular to the x-y plane and at 45° with respect to x-axis (¢ =45°). For the 


random case, the mean beampattern is obtained after averaging the beampatterns for 50 


repetitions of randomly generated array topologies. 











y-axis ideal 
" O position 
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Ww real 
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> 
? A x-axis 
d=N/2 
Figure 7. Ideal and approximately linear arrays with randomly deployed elements and 


position errors. 
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Figure 8. Beampattern of a 10x1 linear array with position errors, compared to a 10x1 


array with equally spaced isotropic elements. Polar angle 6, = 30°, azimuth angle 
g, =45° and mean deviation 30% of 2/2 in both x-y directions. 


Next, the relationship between the position errors and the increase in the average 
sidelobe levels is plotted in Figure 9. The deviation from the ideal position is modeled as 
a uniform random variable in the range of 0 to 2/2 in both x and y directions, L.e., 
mean deviation is 0.54/2. For each value of mean position error, the increase for the 
first sidelobe (red line) and the largest of all the sidelobes (blue line) is calculated by 
averaging 50 simulation runs. It is obvious that the sidelobe power gain increased 
significantly as the mean deviation increased; for example, if the mean deviation is about 
0.44 /2, then the maximum sidelobe increase is almost 5.3 dB while the first sidelobe 
increase is about 4.2 dB. Thus, the strongest sidelobe is lower from the main lobe by 7.7 


dB only while the first sidelobe differs from the main lobe by 8.8 dB. 


19 

















Sidelobe increase (dB) 














L 
10 15 20 25 30 35 40 45 50 
Mean percentage deviation from ideal 4/2 spacing 


Figure 9. Effect of position errors in a random array. Average sidelobe level (dB) for 
the first (red line) and the largest (blue line) sidelobe as a function of the mean 
deviation of the actual element positions from the ideal positions in linear10x1 
array with equally spaced isotropic elements. 


c; BEAMFORMER IMPLEMENTATIONS 


In the previous sections, it was shown that the beampattern of an antenna array is 
determined by the direction of the incoming signal, the number of the array 
elements M x N , and the array topology, which includes the position errors and the set of 


weights w 


ne Lhe objective of a beamformer is to preferentially receive a signal from a 
specific direction or to preferentially transmit a signal in that direction. It is also usually 
desirable to suppress interference signals, which come from other directions. Therefore, 


the beamforming operation consists of adjusting the weights w 


mn 


in such way that the 


main lobe is steered towards the desired signal’s AOA. Such a beamformer for M 
elements in a linear array is depicted in Figure 10, and it is typically used for narrowband 


signals. 
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Figure 10. A M x1 narrow-band beamformer. 


The output of the beamformer is given by 


y(t) = w" x(t) (16) 


where w is the weight vector and x(t) is the signal vector 


xO=[4() 4M .. 4, (OF. (17) 


A conventional beamformer can be described [8] as a delay-and-sum beamformer 
with all weights having the same magnitude. The phases are selected to steer the array 
beampattern towards a desired direction (6,,¢,). However, there are many types of 
beamformers, which can be classified as either data independent or statistically optimum, 


depending on how the weights are chosen [11]. 


In a data independent beamformer, the weights are chosen so as to create a 
specified desired response for all signal and interference cases. The array data (the signal 


vector x(f)) are either not known or not taken into account for the beamforming design. 
If the desired response is F',(@,@) in order to receive a signal from a certain direction and 
cancel out interferences from other directions, then the weight vector w is chosen in such 


a way that the actual response F(@,¢) approximates the desired one. In the following 
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chapters, the beamforming operation will be based on data independent techniques, 


which try to create an approximation of the desired response F,(0,¢) 


The second class of beamformers [11] contains the statistically optimum ones. In 
this case, the weights are chosen based on the statistics of the array data, and the goal is 
to optimize the array response using several criteria. Generally, a statistically optimum 
beamformer tries to cancel out the interfering signals by placing nulls in their incoming 
directions in order to maximize the Signal to Noise Ratio (SNR) at the output of the 
beamformer. Such a beamformer is the Multiple Sidelobe Canceller (MSC), which needs 
auxiliary channels free of the desired signal but is very simple in its implementation. 
Another optimum beamformer requires knowledge of both the desired signal and noise 


covariance matrices R, and R,, respectively, and has the advantage that it maximizes the 


SNR, but it is not easy to implement. The Minimum Variance Distortionless Response 
(MVDR) beamformer computes the weights given a set of constraints and provides a 


very good performance, but it is also computationally intensive. 


There are also adaptive algorithms for beamforming, which compensate for the 
fact that the signal’s statistics are usually not known or may vary with time. Such 
beamformers for the weight determination are based on the well-defined and popular 
LMS and RLS algorithms and also on numerous variations of them. A more thorough 
analysis of these adaptive algorithms and of the previously described statistically 


optimum beamformers can be found in, [11], [12] and [14]. 


D. WIRELESS SENSOR NETWORK ARRAYS 


This section presents a discussion of sensor node deployment schemes and the 
concepts of the processing and communication costs which will be used for the 


performance evaluation of the beamforming algorithms. 


1. Deployment of Sensor Nodes 


As mentioned in Chapter I, recent developments in the MEMS technology have 


enabled the construction of small, cheap, multifunctional sensors with signal processing 
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and communication capabilities. These sensor nodes can be deployed over an area of 
interest and can be used in order to collect process and transmit information. However, 
there are still many challenges for the efficient implementation of a sensor network, such 
as efficient power consumption, controllable deployment of the nodes, source 


localization, self organization, and others. 


In this work, the specific application scenario which is examined is summarized 
in Figure 11. A set of sensor nodes are deployed in the battlefield, and they acquire the 
desired information, which may be any type of signal (acoustic, video, etc.). The UAV, 
flying over the sensor field, establishes a connection with the network in order to obtain 
the collected data. Due to the limited sensor capabilities, a single node can not transmit 
its data directly to the UAV since it does not have the required transmission range. 
However, the nodes can be organized into a large antenna array where each sensor plays 


the role of an antenna element and implements a distributed beamformer. 


Z-aXxiS UAV_-— 








Sensor 
network 





> 


Figure 11. | A sensor network deployment used for information transfer to a UAV (After 
Ref. [34]). 
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In the previous section, it was shown that the best performance for an antenna 
array is achieved when the elements are located on a rectangular grid with an 
interelement distance equal to 2/2. Nevertheless, it is obvious that this ideal deployment 
can rarely be achieved in a real operational environment. The nodes may be dropped by 
the UAV or deployed by a ground force and consequently the sensor array topology 
cannot be the desired one. This randomness imposes position errors, which effect the 


beampattern of the sensor array as analyzed in the previous section. 


Much work has been done in analyzing the dependence of the random array 
response on the statistical characteristics of the array topology [19], [20], [21],[22]. 
Furthermore, much emphasis has been placed on algorithms that allow a node to define 
its position with respect to a reference node within the sensor network [23]. In the next 
chapter, several beamforming algorithms that use location information will be described. 
Assuming that the coordinates of a node in a local coordinate system are known, there are 
many schemes for successfully mitigating the effect of the position errors and choosing a 
suitable set of weights for the beamformer. However, there are limitations in the 


performance of these algorithms when the position errors are large enough. 


Another practical approach could be to search and locate a suitable subset of 
nodes and then form an antenna array. In [24] and [25], a central node (cluster head) tries 
to find a set of nodes whose topology is close to a uniform linear or planar array 
according to some geometric criteria. Using only the distance between the nodes as 
known information, the proposed algorithm tries to find the optimum subset of sensors 
with minimum mean position error. Figure 12 illustrates this concept, where a set of five 
sensors yielding the best approximation to a 5x1 uniform linear array is determined. 
After finding this optimum set of nodes, several beamforming schemes can be applied in 
order to find the weight coefficients. Therefore, the combination of finding a suitable set 


of nodes with an appropriate beamforming technique can provide sufficient performance. 
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Figure 12. Finding five nodes approximating a 5x1 linear array with equally spaced 
elements (After Ref. [24]). 


Zz: Communication and Computational Cost 


An important issue in the sensor networks is the efficient energy consumption by 
the nodes. Low-power hardware components and low-duty cycle operation techniques 
must be applied in order to achieve the required ultra-low-power operation. Energy is 
consumed while processing or transmitting data. Therefore, the beamforming algorithms 
must be evaluated in terms of realistic power consumption. In order to achieve this goal, 
it is obvious that there is a need for schemes implemented in a distributed manner; thus, 
the computational effort and consequently the energy needed is shared among the 
sensors. Moreover, these techniques should minimize the communication since wireless 


transmission consumes considerable amount of power during a node’s operation. 


The need for distributed algorithms in order to minimize the communication 
power consumption is discussed in [26]. A radio transmitting 1 kb of data over a distance 
of 100 m, with an operating frequency of 1 GHz using BPSK modulation having an error 


probability of 10° and fourth-power distance loss with Rayleigh fading, requires 
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approximately 3 Joules of energy. The same amount of energy can perform 300 million 
instructions for a 100 MIPS/watt general processor [26]. This results in a ratio of 30,000 
processing instructions per transmitted bit with equal energy consumption. Other 
practical implementations [26], [27] have yielded ratios from 200 to 3000. This ratio of 
communication cost to computational cost depends largely on the sensor characteristics 
(transmission range, complexity of the instruction in number of bits, etc.) and can be 


increased in the presence of noise where retransmissions will be needed. 


In general, the relationship of the communication and computational cost with the 
power consumption depends on the technical specifications of the sensors, the 
applications, and other factors that are difficult to exactly predict, such as the presence of 


noise. However, a general framework can be formed using the definitions in Table 1. 


a Mean power per transmitted bit (power consumed to transmit a number of 


bits divided by this number) 





Ni Number of transmitted bits 





P Mean power per instruction in the sensor’s processor (power consumed to 


perform a number of instructions, divided by this number) 




















N, Number of instructions 
P. Total transmission power or communication load 
A Total processing power or computational load 
P Total power for a specific application 
Table 1. | Quantities used for defining the communication and computational cost. 
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The transmission power P. is given by 
P=N,,xP, (18) 
and the processing power is given by 

PL,=N,xP, (19) 
while the total power is the sum of these two factors 

P=P.+P.. (20) 


The parameter 77,, is defined as the ratio of the power per transmitted bit to the power per 


instruction, and as mentioned before, it can vary from 200 to 3000, indicating that one 
transmitted bit requires much more energy than one performed instruction by the sensor’s 


processor. 


Ng = (21) 


ro |g 


Throughout this work the communication cost is focused on the implementation 
of the algorithms. Therefore, it depends only on the data elements that need to be 
transmitted for the implementation of the different beamforming schemes. These data 
elements must be encapsulated in data packets, thus there is a communication overhead 
due to these packets. Furthermore, the organization of the sensor nodes in a cluster needs 
also a number of control packets. This networking cost, which depends on various factors 
as the number of sensors, the noise and the protocols used, will not be taken into account 


throughout this work and will be left for future analysis. 


Using the above definitions, the total power consumption for the implementation 
can be calculated and the different beamforming methods can be compared using a 


common framework. 
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E. SUMMARY 


In this chapter, the basics of beamforming in antenna arrays were discussed; 
specifically, the beampatterns of the uniform linear and planar arrays were presented. 
This was followed by an analysis of the effect of the position errors in random arrays on 
the array performance as measured by an increase in the sidelobe levels. Beamformer 
implementations were discussed, including a brief reference to various techniques that 
have been reported in the literature. The specific operational scenario for communication 
between sensor networks and UAV _ was also briefly presented. Finally, the 
communication and computational cost as functions of the power consumption were 
introduced as metrics for the evaluation of beamforming algorithms in sensor networks. 
The strict requirements for low power consumption by the sensor nodes create the need 
for distributed beamformimg algorithms (presented in Chapter IV) compared to 


centralized algorithms (presented in Chapter III). 
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Hl. CENTRALIZED IMPLEMENTATION OF A BEAMFORMER 


This chapter is focused on centralized implementation of the beamforming 
operation for linear arrays. In the centralized approach, it is assumed that all information 
needed to determine the weight vector is available in a specific node, which can be the 
cluster head in a sensor network. This node will collect all the data, such as the steering 
vector d(@,¢), which depends on the array topology, the direction of the desired signal 
(@,,@) and the direction of potential interferences, and will calculate the weight vector 
to provide the desired array response at the output of the beamformer. 

Two different implementations are presented and their advantages and 
disadvantages, such as simplicity in implementation, overall performance and 
computational demands, are discussed. Finally, an analysis of the computational and 
communication cost examines the feasibility of these implementations in sensor networks 


where power efficiency is a major issue. 


A. PHASE MATCH OF THE STEERING VECTOR 


This is the simplest implementation of a beamformer and can be considered an 
extension of the conventional beamformer for the case of random position errors. In this 


method, as in the case of a uniform array (see Chapter II), the weights of the elements w 
are chosen with the goal to match the steering vector d(@,) and create the main lobe with 


maximum gain towards the AOA @,. If the steering vector is 
d(9,) = fl eh sin @% elhs sin i el Pxn sind i (22) 


then the weights are w=d(0,). Each element of the weight vector has unit magnitude 


and the same phase as the corresponding element of the steering vector. Thus, the array 


response 
F(0)=w" d(@) (23) 
has its maximum F(@,) at the AOA 6,. 
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Due to position errors in the sensor array, the distances x,,x,...,x, from the 


reference node are not multiples of the ideal distance 4/2. However, the weight vector 


can still be selected to match the steering vector in the desired direction of 0,. The 


implementation of this approach in a sensor network is very simple and assumes the 


following: 
a) Each node can calculate its position from a reference node accurately. 


b) The relative positions of the nodes are disseminated to a central processing 


node, a cluster head if a cluster hierarchical architecture is established. 


c) The desired steering direction of the incoming signal is known to the central 


node. 


d) Errors due to noise during the communication among the nodes are not taken 


into account. 


The cluster head (central node) collects all the necessary information and 
calculates the weight vector. Then it sends to each node its corresponding weight, which 


will be used for the beamforming and the steering of the sensor array. 


The array pattern of this simple beamformer is shown in Figure 13 for a 10x1 
linear array with uniform position errors up to 0.44/2. The signal’s direction of arrival is 


known (6, =30°, @ =45°). The figure shows the mean beampattern averaged over 50 


simulation runs. 


Some conclusions can be drawn for this beamforming technique. First, it is simple 
in implementation since the central processing node needs only to receive the relative 
position from each node one at a time and then send back the calculated weights to each 
node. Second, the main lobe of the beampattern has the maximum value in the direction 


of arrival (8),@), and it is exactly equal to the maximum value of the (ideal) uniform 


array. However, the beampattern presents notable deviations for angles other than the 
AOA. In this specific case, the sidelobe power level increases significantly for angles 
lower than 0°. In general, the performance of the array approaches that of the uniform 


array at angles around the AOA but deteriorates for other elevation angles. 
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Figure 13. Mean beampattern of a random 10x1 uniform linear array with position 
errors, averaged over 50 simulation runs, compared to a uniform 10x1 array of 
equally spaced elements. The AOA is (6, =30°, ¢ =45°). The position errors 


are uniformly distributed between 0 and and 0.44/2 in both x and y directions. 
Beamforming is implemented by matching the steering vector. 


B. BEAMPATTERN APPROXIMATION IN THE LEAST SQUARE SENSE 


The beamformer presented in the previous section computes the weights of the 
antenna elements by trying to achieve the best performance in the direction of the signal 
arrival and does not set any requirement for the other directions. Although it behaves well 
at angles around the desired direction, it suffers from high sidelobes in other directions. 
Therefore, it cannot be used when there are interference signals coming from directions 
other than the signal’s AOA or if the desired signal comes from a certain range of 
directions. If there is a strong source of interference, the desired array response should be 
zero in that direction, in order to cancel the interfering signal. As a result, there is a need 


for a beamforming design, which calculates the weight coefficients w in such way that 


the resulting response approximates the desired response over a range of directions. 
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1. Least Squares Problem Formulation 
As mentioned before, the array response for a specific angle 9 is given by 
F(@)=w' d(). (24) 


If the desired response F',(@) is defined over m number of angles, then the actual array 


response F(@) of n antenna elements for these m angles is 
F(0)" =w" D6) (25) 
where 
FQ@)=[F@) FO) ~- FO)) (26) 
is a mx1 vector which contains the array response for each one of the angles @,, with 


i=l,...m. The weight vector w is an nx|column vector and D(@) is an nxm steering 


matrix containing the steering vectors d(@,) for each one of the angles @,: 


D(O)=[d(@) dA) -- d(4,)] (27) 
eh sin), eh sin 0, Ne! elPx sind, 
JBx sinh eh sin 0, Sos el Px sind, 
or DO)=| ; ; (28) 
el Pxn sin}, elPxn sin 0, wee e/Pxn sin8,, 


Note that for simplicity the array response is defined over a set of arrival angles 


6,,...,8 


m 


while the azimuth angle ¢, is fixed. If the desired response is defined over a 
combination of m elevation angles 6,,...,0, and g azimuth angles Dy s++9Pys then the 
array response F(6,¢) is a mxq matrix and the steering matrix D(@,¢) isan nxmxq 


3-D matrix. This adds unnecessary complexity to the formulation of the problem and is to 
be avoided. In the following analysis, the array response will be defined for a known, 


fixed azimuth angle ¢, and a set of m arrival angles 6,,...,0,,. 
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If the desired response is F'_,(@), the goal is to select the weights w such that the 
actual response F(@) generated by this set of weights approximates the desired response. 
This can be done by trying to minimize the L,-norm of the difference between the 
desired and the actual response. In other words, the weight coefficients w are chosen in 


order to minimize the mean squared error between F',(@) and F(@) 
e =min|F(0)-F,(0) . (29) 


Therefore the problem of the antenna weight calculation is solved using a classic least 


squares (LS) procedure. Taking the transpose of (25), the array response can be written as 
F(0)= D@)" w. (30) 


For m number of angles or approximation points and n number of sensors, with m>n, 


the problem becomes an overdetermined LS problem 
e=min|D(0)" w- F,(0)| 31) 


where again F(O) isa mx1 vector and D(@)" is a mxn matrix. Provided that the nxn 
matrix D(@)D(@)" is invertible (i.e., D(@)" is full rank), then the solution to the LS 


problem of (30) is given by 

w=D'(@)F,() (32) 
where D*(@) is the pseudo inverse of D(@)" defined as 

D' (0) =(D(@)D()") 'D@). (33) 


The weight coefficients that have been calculated using the LS approach have both 
different amplitudes and phases. In the array with equally spaced elements, the 


amplitude /,, was the same for all elements (/,, = 1,Vm) and only the phase was different. 


However, in the LS approach the weights have also different amplitudes. So in the LS 
approaches the sensors modify both amplitudes and phases in order to approximate the 


desired response. 
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The design of an array using the above decribed LS approximation of the desired 
array response is related to the design of arrays using the Fourier analysis. In this 
approach, the desired array function, which is required to create the desired array 
response, can be expanded in a Fourier series. Truncating this Fourier series, the result is 
an array with a finite number of elements. Such an array is optimum in the sense that no 
other array with the same number of elements can approximate the desired array function 


with lower mean squared error. 


2. Implementation and Performance Analysis 


The significant advantage of the above LS approach is that the desired response 
can be closely approximated over a set of angles, thus providing the ability to cancel out 
unwanted signals and amplifying only the desired one. Results of such an implementation 
of the LS scheme for a 10x1 random array can be seen in Figure 14, where the desired 
response F',(@) is the array response of the (ideal) uniform 10x1 linear array. For the 
approximation of the desired beampattern 180 points (angles from —90° to 90°) were 
used. The mean beampattern is averaged over 50 simulation runs, and there is virtually no 
difference between the desired (red line) and the actual beampattern (black line) as 
generated by the weights computed by minimizing the LS criterion. The actual response 
can be also compared with the response derived by the technique of the phase match 


approach from the previous section (cyan line). 
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Figure 14. | Mean beampattern of a random 101 linear array with position errors (blue) 
and mean beampattern calculated using the “phase match” (cyan) and the LS 
approach (red), averaged over 50 simulation runs, compared to a uniform 10x1 
array (black). The AOA is (8, =30°, @ =45°). The position errors are uniformly 


distributed between 0 and and 0.52/2 in both x and y directions 


The performance is almost the same for both phase match and LS schemes around 
the desired AOA, but the LS approach is obviously better in any other direction. The 
phase match technique is more vulnerable to interference signals, especially to those 
coming from @<0°, while the LS approach creates a beampattern almost identical to the 
desired response of a uniform array. Indeed, the determination of the weights using the 
LS approximation technique achieves elimination of the effect of the position errors (blue 
line); the position errors are uniformly distributed between 0 and 0.54/2 in both x 


and y directions. 


The LS approach is clearly a centralized approach. The implementation in a 
sensor network makes similar assumptions as in the phase match approach including an 


additional assumption: the desired array response defined over a set of angles 6,,...,0,, is 


m 


33 


assumed to be known to the cluster head. This response may include the direction (@),4,) 


of the incoming signal and/or the directions of interfering sources. 


The cluster head collects all the necessary information and specifically the 


relative positions x,,...,x, of the nodes and the desired response F',(@). The steering 
matrix D(@) is completely defined by the positions and the desired angles for the 


beampattern approximation. After collecting all necessary data, the cluster head 
constructs the steering matrix using the position data and the set of angles and calculates 
the weight vector by solving the LS problem of (30). Each node receives the 
corresponding coefficient from the cluster head, and the sensor array is ready to transmit 
or receive towards the desired signal’s direction and simultaneously suppress potential 


interfering sources. 


The LS approach’s performance is satisfactory in eliminating the effect of the 
position errors as observed in Figure 14. The assumptions mentioned before, such as that 
information about the relative position of the nodes and the desired response (incoming 
direction of signal and interferences), are reasonable in a sensor network environment. 
There are algorithms which allow a sensor node to calculate its relative position from a 
cluster head using localization techniques [23], [26]. Thus, the LS scheme could be 
implemented in a cluster-based sensor array where all information is collected and 


processed in the cluster head. 


Nevertheless, this centralized approach suffers from the inherent problems of any 
scheme where the processing effort is not distributed among the nodes but performed by 
a single node (cluster head). Since the main effort is undertaken by one sensor only, a 
single point of failure is created in the sensor array. If this specific node fails, then the LS 
problem has to be solved from the beginning, i.e. a new central sensor must collect all the 
information, calculate the weight coefficients and disseminate them to the other nodes. 
This is a waste of valuable processing and communication resources which are very 
restricted in the sensor nodes. Considering that failure of sensor elements in a sensor 
array is something common due to their low battery life and low-cost construction, this 


centralized approach lacks robustness. 
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Furthermore, since new nodes may offer to participate in the sensor array or other 
nodes may decide to switch into sleep mode due to limited power, the array topology will 
change frequently and dramatically. This means that, for any modification in the sensor 
array, the LS problem must be solved from the start in order to determine the new set of 
weight coefficients. For this reason, the centralized approach does not offer the required 


scalability, which is desired in a sensor network environment. 


Finally, the processing load of solving the LS problem increases exponentially as 
the number of sensors n and the number of approximation points (angles) m increase. If 
this processing effort is performed by a single sensor, then its power will be consumed 


rapidly and the sensor will fail. 


The solution is a distributed implementation of the LS approach, and such two 
schemes will be presented in Chapter IV. First, the single point of failure limitation can 
be overcome by distributing the processing load among many nodes. The total processing 
effort in a distributed approach may be higher compared to the centralized approach, but 
the tradeoff is the increased robustness of the system. Furthermore, the LS approach must 
be implemented in such way that changes in the array topology do not require the 
solution of the problem from scratch but only small modifications to the already 
calculated weight coefficients. This can be achieved by an algorithm implemented in a 
distributed fashion, which will offer considerable scalability to the problem. Since the LS 
approach performs very well in approximating the desired pattern, a distributed 
implementation of it will first retain this feature and second, it will offer the desired 
scalability and robustness; thus, it can provide a reliable solution for beamforming in 


sensor networks. 


In the next section, a discussion about the processing cost as a function of the 
number of approximation points is presented. As mentioned before, the processing load 
increases dramatically with the number of approximation points; therefore, the minimum 


number of points should be used in order to reduce the computational effort. 
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C. APPLICATION OF LEAST SQUARES BEAMFORMING APPROACH TO 
SENSOR NETWORKS 


In the previous section the LS approach was presented for the approximation of a 
desired response defined over a set of direction angles or approximation points. As 
mentioned before, the processing load of the LS approach depends largely on the number 
of nodes and the number of approximation points. For example, if 1 is the number of 


sensors and m the number of angles, the solution of the LS problem using the normal 
equations needs (m+"% )n? flops [28]. By using fewer approximation points, the 


processing load can be reduced, which is crucial for the sensor nodes with limited 
processing power. However, the solution will not be the same as the unmodified case, 
and the approximation will degrade as the number of points is decreased. Therefore, it is 
worth examining the “quality” of the desired response approximation as a function of the 


number of approximation points. 


In order to illustrate this relationship, a series of 50 simulations with 
corresponding random arrays was performed. In each run, the desired response of a 
uniform linear 10x1 array (ideal) was approximated using different number of 
approximation points ranging from 10 to 90. For each set of approximation points, the 
corresponding solution for the weight coefficients was derived using the LS approach. 
The reference approximated response and the reference weight vector was based on 180 
approximation points (one for each degree from —90° to 90°). The solution of the LS 
approach for each set of approximation points is compared to the reference solution 


derived from 180 points. 


Two error metrics are used to evaluate the performance of each set of 


approximation points. The first metric e,, is the mean L,-norm of the difference between 


the reference weights and the examined weights 


138 
es W;—Wigo| (34) 





38 


where w,,, are the weights derived from the 180-point approximation (reference), w, are 
the weights derived by an i- point approximation, for i=10,...,90, and S =50is the 
number of simulation runs. 

The second metric e, is the mean L,-norm of the difference between the 
reference array response derived from 180 approximation points and the array response 
for fewer approximation points (10 to 90) 

i 


er = TYR Feo (35) 


where F,,,(@) is the approximated array response using 180 points (reference) and F,(@) 
is the response derived by an i- point approximation, for i=10,...,90. 
The test topology is a 10x1 random array where the antenna elements have 


position errors, modeled as uniform random variables between 0 and 0.4/4/2. The results 


are shown in Figures 15 and 16 for the two errors e, and e,, respectively. 
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Figure 15. Error metric e, of weights as a function of the number of approximation 
points m. 
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Figure 16. Error metric e, of approximated array responses as a function of the number 
of approximation points m . 


In Figure 15, error metric e,, decreases as the number of approximation points 


increases, which is reasonable. The significant result comes from the fact that the 
difference between the reference weights and the weights computed for fewer 
approximation points is very small, even for 15 approximation points. This means that 
using only 15 approximation points, the weight vector is very close to the reference 
weight vector (180 approximation points). Therefore, there is no need to use 180 points to 
approximate the desired response since this can be done by using only 15 points, which 


implies that the processing cost can be significantly reduced. 


Similar comments can be made about the results in Figure 16 for the error 
between the reference array response derived by 180 approximation points and the 
responses derived by fewer points. If the desired response is known and needs to be 
approximated using a LS approach for the determination of the weights, the number of 


approximation points can be substantially lowered without significant degradation in the 
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performance of the antenna array. Thus, the processing cost can be reduced, which is 


desirable in a sensor network environment with restrictions on the power consumption. 


D. SUMMARY 


In this chapter, two centralized approaches for the computation of the weights 
were presented. The first technique is simple and easy to implement, but the performance 
is not satisfactory in the presence of interference signals. The second approach is based 
on the LS approximation of a known desired response by selecting the weights in order to 
minimize the error between the actual and the desired response. Although having 
satisfactory performance, the LS scheme still lacks robustness and scalability since it is a 
centralized approach. Finally, the performance of the approximation of the desired 
response as a function of the number of approximation points was examined. The next 


chapter presents two algorithms that implement the LS approach in a distributed fashion. 
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IV. DISTRIBUTED ALGORITHMS FOR BEAMFORMING 


A Least-Squares (LS) based centralized approach for solving the beamforming 
problem was described in the previous chapter. Nevertheless, in a sensor network 
environment where the processing and communication load must be shared among the 


nodes, there is a need for distributed beamforming algorithms. 


Two distributed implementations of the LS beamforming problem are proposed in 
this chapter. The first is based on the well known QR factorization using Householder 
transformations, where each node performs a part of the QR decomposition process. The 
second approach is a distributed iterative solution of the LS problem, which converges 
quickly to the actual weight coefficients. Both algorithms efficiently distribute the 
processing load among the nodes; however, the tradeoff consists of an increase in the 


communication cost. 


A. DISTRIBUTED QR FACTORIZATION WITH HOUSEHOLDER 
TRANSFORMATIONS 


1. Householder Transformations 


A Householder transformation (or Householder reflection) is a transformation of 


reflecting a vector about some known plane [28], [29]. Given an arbitrary vector 








m 














ed pe He. se a | and a unit vector e, =I 0 -: 0] < Pix © thie 








mxm Householder matrix is defined as the matrix that transforms x to a vector parallel 


to é, 
Hx=ae, (36) 
where a is ascalar. 


A Householder matrix can be formed by [29] 


H =I" (37) 
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where f is a scalar, J is an identity matrix and v is an mx1 vector. From the above 
definition, the Householder matrix H is completely determined by vector v and scalar 
£, so there is no need to store all m* elements of H, but only vy and, which are only 


m+1 elements. 


The significant characteristic of the Householder matrix is that it is an orthogonal 


matrix; thus, it has the property 
H'=H'=H (38) 


Orthogonal matrices or orthogonal transformations can be used to obtain a QR 
factorization of a matrix A, and this in turn can be used to solve a linear system Ax =b, 


as described in the following section. 


The above defined Householder matrix zeros out all the elements of an mxl1 


vector x except the first one, but one can construct Householder matrices that zero out 
only the last m—k components of a vector x [29]. Let x) and x?) define 


1 
x 


=[% Ny ale (39) 


x =[x, a os Ail (40) 


while I and 1 denote (k-1)x(k—-1) and (m—k+1)x(m—k +1) identity matrices, 
respectively. From (32), an (m—k+1)x(m—k+1)Householder matrix H{” can be 


constructed as [29] 
He =1 -—_y,v, (41) 


such that 


Hx? =|2| e, . (42) 


2 


By defining the mxm matrix H, as 
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H, a : | (43) 


the following results 


T 
JO x x mn 1/2 

Hx Sl = | apeaa pe oe Di 0. O|. (44) 
fe ago Ay x. i=k 


The Householder matrix H, defined in (43) is an orthogonal matrix that can be 


written as 


1 
A, ieee (45) 


k 





m 





and 





and acts like an identity matrix on the first k—1 coordinates of any vector x € 





transforms the rest into a unit vector. Moreover, it is not necessary to store the entire 


matrix H, , but it is enough to store the (m—k+1)x1 vector v, and the scalar 7, . 


2. QR Decomposition 


The QR decomposition is used in order to factor a matrix A into a product of two 


matricesQ and R [28], [29] 
A=QR . (46) 


If A isan mxn matrix, then Q is an mxm orthogonal matrix and R is an mxn upper 
triangular matrix. The QR decomposition is used to solve a linear system Ax =), and it 


can be computed by applying a series of Householder transformations to matrix A. 


By defining 
A, A, a3 a, 
Gy, yy Ag a>, 
A=| 43, Gy G3 a3, (47) 
| ni ang ang Qin 
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and using the procedure described in the previous section, an mxm Householder matrix 


H, can be found, which when applied to the first column of A will give a multiple of e,. 


Thus, the result of multiplying H, with A will be 


re) (1) (1) (1) ] 
GA, A Az co a, 
() (1) (1) 
0 ay ay * G, 
as () (1) (1) 
MA=| 0 ay a3 ++ G;, (48) 
(1) qd) qd) 
0 ano an3 Ginn | 








where the superscript denotes that the matrix element is affected by the 4H, 
transformation. Note that H, is fully described by vector v, and scalar £, as mentioned 
before. Then, a second Householder transformation H, that will zero out the last n—2 


elements in the second column of H,A can be again constructed. Thus, H,H,A will be 


of the form 
[4d qd) qd) (1) 
GA, Ay 43 7 GA, 
(2) (2) (2) 
O° sa ay Fo, 
H,HA=|0 0 a? -- a® (49) 
0 0 a of] 








where again the superscript indicates the last transformation that affected the 
corresponding element. For example, the first column and the first row are not affected 


by the H, transformation. From (41), H, is of the form 


I, 0 
se 0 H® (50) 
where J, is a 1x1 identity matrix, and H{” is an (m—1)x(m-—1) Householder matrix. 


Since H{” is fully determined by an (m—1)x1 vector ve and a scalar £,, only 


m elements are needed to be stored for H”’. 
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Continuing to apply the Householder transformation in this way will result in an 


upper triangular matrix R, ie., 








H,H ,--H,H,A=R (51) 
where R is of the form 
arias ae ai, 
0 a af a? 
0 a ae 
0 
R= (52) 
0 
10 0 0 0 Of 


The mxn matrix R is composed of two matrices, a nxn upper triangular R, and a 


(m—n)xn zero matrix: 


R, 
r-|5 | (53) 


OF =H tht, =| 4 (54) 


where Q, is an nxm matrix containing the first n rows of Q’, and Q> is an 
(m—n)xm matrix containing the last m—n rows. Since Q'A=R and Q is orthogonal, 
it follows that 


R 
A=OR=[Q, 2.) |=: (55) 
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The above derived QR decomposition can be used in order to find the solution to 


the least squares problem Ax =b, where A is an mxn matrix and x and b are nx1 and 


mx1 vectors, respectively. The solution x satisfies the normal equations [29]: 

A’ Ax=A'b (56) 
and by using the QR decomposition, they can be written in the form 

R'QOARx=R' Ob . (57) 


By setting 


rp, [QS] _|S 
-oe-[2).(5 ie 


and taking into account that QQ, =I and R_ is nonsingular, (57) simplifies to 


R=, (59) 
The system of (59) can be solved using back substitution [29] and the solution 

= Ri'c, (60) 
is the unique solution to the least squares problem. 


In summary, if A is an mxn matrix, the least squares problem can be solved 


using the QR decomposition by Householder factorizations as follows: 


a) Compute R=H,H, ,:--H,H,A and c=H,H,,::-H,Ab. 


R é 
b) Partition R and c into a block form R -| | and c 2 where R, and c, 
Co 


each have n rows. 


c) Finally, solve R,x =c, using back substitution. 
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3. Proposed Algorithm Description 


In Chapter III, a method for determining the weights that provide the best 
approximation of a desired response in the LS sense was described. The LS problem of 
(30) can be solved using the QR factorization presented in the previous section. The 


cluster head collects all the positions from the sensor nodes and constructs the steering 
matrix D(0)" . Then it computes the QR factorization using Householder transformations 
and finds the solution of the weight vector w for a given desired response F,(@). The 


procedure is straightforward, but it lacks robustness since a single node bears the entire 
processing load, which increases the possibility of failure, causing the problem to be 
solved again from scratch. Therefore, a distributed algorithm for solving the LS problem 
with the QR decomposition is desirable. The processing load is shared by the nodes, but 


there is also a tradeoff of increased communication load. 


For the implementation of the algorithm, the following reasonable assumptions 


are made: 

a) Each node can calculate its position from a reference node accurately. 

b) Each node can broadcast information to other nodes in the cluster. 

c) The desired array response for a set of m directions @,...,0,, of the incoming 
signal is known. 


d) Errors due to the noise during the communication among the nodes are not 


taken into account. 


The algorithm exploits the specific nature of the matrix D(@)" to be decomposed 


by the QR factorization. By taking the transpose of (28), the steering matrix is 


j Px, sin O, j Bx sin B, j Px, sin), 
elf 1 1 el 2 1 one el n 1 
J Bx, sin 6, e/P» sin 0, a3 eiPxn sin 8, 
H 
DO) =|": (61) 
IPN SinOn — IPSN |, p JB, SiNOy 
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where n is the number of sensor nodes and m the number of angles or approximation 


points. Obviously, the first column is constructed by the position of the first node x, and 


the set of angles 6,,...,0,,. Similarly, the elements of the second column depend on the 


second node position x, and the set of m angles. In general, each column i of D” (0) 
depends only on the position x, of the corresponding node and the desired angles. 

From (61), it is obvious that since one node knows its position from a reference 
node and the set of m angles, it can construct its corresponding column of D” (@) 


without exchanging any information, such as their position, with the other nodes. Thus, 


the matrix D”(@) can be stored in a distributed way among the sensor nodes. 

From the analysis of the previous section comes the observation that the first 
Householder transformation H, is needed to zero out all the elements of the first column 
except for the first element. Since for the construction of matrix H, only the data values 


from the first column are used, the conclusion drawn is that it can be computed by the 
first node only. Therefore, the first node, which initializes the QR decomposition, 


calculates the H, matrix and multiplies its own first column by this matrix. Then it 
broadcasts H, to all other nodes, which in turn transform their corresponding columns 
with H,. There is no need to transmit the entire matrix H, but only the vector v, and the 
scalar f£,. After this first step of H, computation, transmission and application is 
completed, the distributed steering matrix will have the form of (48): 


[ () (1) (1) (1) 
ay, 


ay Oia ie Ae 
(1) (1) (1) 
O ay ay + A, 
H = H = 1 1 1 
Di (9) = HD" (@)=| 0 as ay as, (62) 
(1) (1) (1) 
0 and an3 _ Gann | 








ad 


where ay are the matrix elements transformed by H,. 


The second step is similar to the first and is performed by the second node. The 


second node computes the H, transformation, which depends only on the second 
50 


column, and then H, (vector v, and scalar £,) is broadcast to each node which in turn 


applies the new Householder transformation to its column. The resulting distributed 


matrix is 
[ (1) (1) (1) (1) | 
Ay Ay Az " GQ, 
(2) (2) (2) 
O ay Gy + Gy, 
H = H _ 2 2 
D> (8) = H,H,D (A)=| 0 0 Ge a a? . (63) 
0 0 ak a2 








From the structure of D5(8), the following observations are made. The first 


node receives the H, transformation, but it does not need to apply it to its column since 
H, does not affect the first column. Additionally, each node does not need to apply the 
H, to the first element of its column since the first row is not affected by the H, 
transformation. Similarly, the following Householder transformations are applied only 


when needed, thus avoiding redundant calculations. 


The procedure goes on until all Householder transformations 4H,,...,H, are 


computed, broadcast and applied by the nodes. The distributed matrix will result in an 


upper triangular form 


ae) (1) () (1) 

ay Ay 43 a, 

(2) (2) (2) 

0 ay a5; a5, 

(3) (3) 

0 ay a3, 
D" (0)=H.H H,H,.D" (6)= ‘ : 64 
‘ny ( )= ntt nl" yeaa | ( )= . (n) ( ) 

Qin 

0 








10 0 0 0 0] 
The QR decomposition of the D”(@) matrix is achieved in a fully distributed way, and 


the processing load is shared among the nodes. There is no additional processing effort 


since the Householder transformations are applied only to the affected elements, as 


a1 


described before. There is no difference between the centralized and the distributed QR 
decomposition in terms of the processing load. The significant advantage is that instead 
of one node bearing all the computational burden, each node in the array takes its turn in 
completing the QR decomposition. However, there is some additional communication 


load since the H, matrices need to be broadcast. This is the tradeoff for relieving the 


cluster head from the heavy processing load. This first phase of the algorithm is depicted 


in Figure 17. In the 1° step the 1 node calculates the matrix H, and broadcasts it to all 
the other sensors. In the 2™ step the 2™ node calculate the matrix H , and the procedure 


continues until the last node which computes the last matrix H, 


1“ step 


@ 


OF 
\ 
© 


Hi 
1, H2 


oO 02 O 


Hi, H2, H3 





M1, H2 H1, H2 Se 
o) ©) H3 


Figure 17. First phase of the algorithm for distributed QR decomposition by the 
sensor nodes. Bolded and underlined H, above the nodes denote the Householder 


transformations that are already stored in the sensor. Simple H, denote the just 
computed and broadcast Householder transformation. 


The next phase in solving the LS problem is to update the desired response using 


a series of Householder transformations to obtain vector c: 
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: FG) © 
F;”(0;) 
F;(6) 
=H Ag HAF O)=|) 2 (65) 
F;"(6,) 
LF) ,.) | 
C 
which is then partitioned to c =| — | according to (58). 
we) 
The third phase includes the solution of the system 
Cc 
C 
KWH CS) (66) 
c 


by back substitution where R, defined by (52) and (53) is an upper triangular nxn 


matrix containing the first n rows of Di (9) ; 


[ (1) (1) (1) () | 


a, Ay Az « GQ, 
(2) (2) (2) 
Oy "Gay 28% Gs, 
. 3 3 
R,=| : 0 as vee a : (67) 
0 0 O a” 








The n" node needs to solve the last equation of the system 
aw =c (68) 
and immediately computes its weight w, for beamforming. The (n-1)” node needs to 
solve the (n—1)" equation in order to find its weight w,_,: 


GW a Fa WSC ae (69) 


n—-l,n-1"" n—- n-l,n'"n 
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The a‘";? , data element of matrix R, in (67) is constructed by the (n—1)" node by 


n-l,n-1 


— pIPXy1 Sin 8,4 
n-l1,n—1 ~ e 


applying all transformations H, toH, , to the initial element a since 


(n-l) 


n-l,n ? 


this data element belongs to its own (n—1)" column. However, it needs the valuea 
which is in the n” column; thus, it depends on the position x, of the n" node. It also 
needs to receive the weight w, and compute the value c, ,. The n™ node broadcasts its 
position x, and its weight w, to all nodes. Now, the (n—1)” node can solve (69) and 


calculate the weight w, ,. 


Similarly, the (n—1)" node broadcasts its position x, , and its weight w,_, to all 
the other nodes, which need this information to solve their own equation. For example, 


the (n—2)" node has to solve the equation 


(n-2) 
n—-2,n- 


(n-2) 
n—-2,n- 


(n-2) aa 
Wy tOno.aWn =Cn-2> (70) 


n-2,n°"n 


as SW ore 


(n-2) 
n—-2,n-1 


(n-2) 
n—-2,n 


so it needs the positions x, and x, _, in order to construct the a and a elements, 
respectively, by applying all the Householder transformations up to H,,,. It also needs 
the weights w,, and w, from the previous nodes. 

The final result is that every node obtains its own weight, and the sensor network 
is now ready to cooperate in order to form an antenna array. The third phase of the 
algorithm is depicted in Figure 18. During the Ist step the last node (3™), calculates its 


own weight w, and broadcasts it along with its position x,. In the 2" step the 2"* node 
uses the received information to solve for its weight w,. Then it broadcasts w, and its 


position x,. The procedure goes on until the 1“ node calculates its own weight w, 
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1“ step 


Q) “2 x3, <3) 


Last phase of the algorithm to implement distributed back substitution by 


the sensor nodes. Bolded and underlined x,,w, above the nodes denote the 
positions and weights that are already stored in the sensor. Simple x,,w, denote 
the broadcast position and just computed weight. 


Computational and Communication Cost Analysis 


The computational cost of the procedure presented in the previous section can be 
measured in number of instructions needed for the implementation. It is known that the 


number of flops for the solution of the LS problem in a central processor is given by [28] 
N, =2n*(m—n/3)+mn+n° (71) 


where the first term is for the determination of the QR factorization, the second term for 
the update of the desired response vector with the Householder transformations, and the 


last term accounts for the back substitution to solve for the weights. 


The significant characteristic of the proposed distributed algorithm is that, as 
described before, no redundant computations are made. The Householder transformations 


are calculated in exactly the same way as the centralized approach, which yields the same 
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number of computations. In the centralized approach, however, each element of the 
distributed stored steering matrix Di)(9) , or similarly the upper triangular matrix R, of 
(67), is computed by a single sensor and no redundant computations are made. For 
example, in the distributed approach, the first node does not construct the entire matrix. 
In the first phase, it uses only the first column for defining the first Householder matrix 
H,, and then in the back substitution phase it calculates only the first row of the R, using 
the H,. Similarly, the second node calculates H, using the second column and updates 


the second row, which is used for the back substitution. The following matrix shows the 


elements of R, with respect to the node that uses and calculates them: 


RSS) OF By eS (72) 








[9 0 0 ++ n{ 
Therefore, no additional computations are made, and the total processing cost is exactly 
equal to that of the centralized approach given in (71). Thus, without any increase, the 


processing load is shared among the nodes, and the cluster head is relieved from the 


heavy computational effort. 


The number of instructions N, can be used to express the processing power using 


the definitions of Chapter II for computational cost. Therefore, from (19) and (71), the 


processing power P, is given by 
P, = N,x P.=[2n°(m—n/3)+mn +n? |x P ; (73) 


However, this distribution of the processing load comes with a tradeoff, which is 
an increase in the communication load. Note that throughout this work the 
communication cost is only the number of data elements that need to be transmitted for 
the implementation of an algorithm. The transmissions required for the coordination of 
the sensors in the distributed algorithm are not taken into account for the calculation of 


the communication cost. In the centralized approach, the cluster head collects all position 
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data from n nodes and sends back n weights, so the number of transmitted data 


elements N, is equal to 


N,=2n. (74) 


t 

On the other hand, in the distributed algorithm, the number of transmitted data 
elements is calculated as follows. The first node sends the H, matrix, or equivalently the 
mx1| vector v, and the scalar £,, a total of m+1 data elements. Similarly, the second 
node broadcasts the matrix H, or the (m—1)x1 vector vy, and the scalar /,, a total of 
m data elements. This procedure is continued until the (n—1)” node, which transmits 
the H, ,. In general, an arbitrary node i transmits the matrix H,, te. (m—itl) 
elements for v, and one for £,. The total number of transmitted elements by all nodes for 


the first phase of the QR decomposition is given by 





ie =(m+2—n/2)(n-1) . (75) 


Ci SS Gein) Greed 


For the phase of the back substitution, the n” node transmits its position and its 
weight, a total of two elements. Similarly, each node except the first one, broadcasts its 
position and weight to the other nodes. Thus, the total number transmitted elements for 


the back substitution is given by 
Cys =2(n—-1) . (76) 


Finally, the total number of data elements to be broadcast during the implementation of 
the algorithm is 

C=Cop t+ Cos =(m+4—-n/2)(n—-]1). (77) 

Using again the definitions from Chapter II for the communication power, the 


power consumption due to communication can be derived. Assuming that each element is 


represented by b bits, the number of transmitted bits is 
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N,,=Cxb (78) 


and from (18) the communication power P. is given by 


P.=N,,xP, =(m+4-n/2)(n-))xP, xb. (79) 


Therefore, the total power for the implementation of the algorithm is the sum of (73) and 


(79) 
P =| 2n?(m—n/3)+mn+n? |x P.+(m+4—n/2)(n-1)x P, xb (80) 


and depends highly on the specific characteristics P,P, andb of the sensors. Similarly, 


the total power for the centralized implementation is calculated to be 
P =| 2n?(m=n/3)+mn+n? |x P +2nxP, xb. (81) 
The following figures summarize the results and compare the two implementations. 


In Figure 19, the number of required instructions N, is plotted as a function of the 
number of approximation points m for a 10x1 array of sensors while in Figure 20, N, is 


plotted as a function of the number of sensors n, for m = 20 approximation points. The 


plots for the processing power P, can be obtained by multiplying the number of 


instructions with P. as in (73). Another important point is that the main processing effort 


happens during the QR decomposition phase as the back substitution procedure is not so 
demanding. This is significant because the QR decomposition is computed once and then 
the results are kept for future modifications. For example, if a new array response is 
required, it is only necessary to compute the updated vector c, of (65) and then solve the 
system of (66) by back substitution; R, does not change. Therefore, in the case of 
communication with a UAV, which is moving and causing the desired response to change 
continuously, the sensor array only has to quickly perform the back substitution phase in 


order to compute the new weights. 
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Figure 19. Processing cost of the distributed algorithm as a function of the number of 


approximation angles and for n=10 sensors. Multiplying N, by P. gives the 
required processing power. 
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Figure 20. Processing cost of the distributed algorithm as a function of the number of 


sensors and for fixed number of approximation angles m= 20. Multiplying N, 


by P. gives the required processing power. 
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The communication cost as a function of the number of approximation points and 
the number of sensors is depicted in Figures 21 and 22, respectively. Figure 21 shows the 
communication cost for a 10x1 sensor array in terms of the transmitted elements. 
Assuming that each data element is represented by b=32 bits in a single precision 
floating point representation, the total transmission power can be found by multiplying 


the number of elements C by P, xb as in (78) and (79). Similarly, in Figure 22, the 


number of transmitted elements C is plotted as a function of the number of sensors n 
for m= 20. Compared to the centralized approach, the communication load is increased, 
especially due to the first phase of the QR decomposition. However, as described before 
for the processing load, this phase has to be implemented only once. Therefore, any 
change in the desired response needs only the back substitution phase, which requires a 


considerably lower number of elements to be transmitted. 
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Figure 21. Communication cost of the distributed algorithm as a function of the 


number of approximation angles and for a fixed number of sensors (n =10). 
Multiplying the number of data elements with P,xb gives the required 


transmission power. 
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Figure 22. Communication cost of the distributed algorithm as a function of the 


number of sensors and for fixed number of approximation angles (m= 20). 
Multiplying the number of data elements with P,xb gives the required 


transmission power. 


In Figures 23 and 24, the total power consumption as a function of the number of 
approximation points and the number of sensors, respectively, is plotted. Assuming that 


P, is equal to one unit of power, and using the ratio 77,,, the factor P,, can be substituted 


by Fi, =7,, xf. Dividing the total power P of (80) by P yields the normalized power 


P. in terms of the required units of power P.: 


P. == [20° m—n/3) + mn +n? ]+(+44=n/ (0-1 xbx 1, (82) 


From the figures, it is obvious that the distributed algorithm requires more power 
than the centralized one, and this is caused by the additional communication load. 
However, this total power cost is shared among the nodes of the sensor network and is 
not undertaken by a single node. For example, Figure 24 shows that a deployment of 


twenty sensors needs about five times more power for the distributed approach than the 
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centralized approach. Therefore, the cluster head in a sensor node consumes four times 
more power in the centralized scheme than the average sensor nodes consume in the 
distributed one. This obviously will cause the cluster head to fail quickly, which means 
that a new cluster head will need to be selected and all the computations from the 
beginning will need to be re-done. Thus, the increase in the power consumption in the 
distributed approach can be considered reasonable if the robustness of the network is a 


crucial requirement. 
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Figure 23. Normalized power P, (number of P.) for both distributed and centralized 


approaches as a function of the number of approximation angles for a fixed 
number of sensors (n = 20). 7,, is assumed to be 200 and b = 32 bits. 
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Figure 24. Normalized power P, (number of P.) for both distributed and centralized 


approaches as a function of the number of sensors for a fixed number of 
approximation angles (m= 20). 77,, is assumed to be 200 and b = 32. 


B. DISTRIBUTED ITERATIVE SCHEME FOR SOLVING THE LEAST 
SQUARE PROBLEM 


In this selection, a distributed scheme based on an iterative solution of the LS 
problem is presented and evaluated. The iterative procedure is performed in steps by all 
the nodes in the sensor array. Starting from an arbitrary initialization for the weights, this 


method quickly converges to the actual solution of the LS problem. 
For the implementation of the algorithm, the following assumptions are made: 
a) Each node can calculate its position from a reference node accurately. 
b) Each node can broadcast information to other nodes in the cluster. 


c) The desired array response for a set of directions @,...,0,, of the incoming 


signal is known. 
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d) Errors due to noise during the communication among the nodes are not taken 


into account. 


The linear system of (31), which is to be solved in the LS sense, is considered a 


minimization problem of the following form: 





€=min Ax—b| (83) 
where 

A=D(0)" 

x=w (84) 

b=F,(8) . 


1. Proposed Algorithm 

The algorithm is based on various parallel methods proposed in the literature for 
solving the LS problem [30], [31], [32], [33]. 

Let A be an mxn matrix while x and b are nx1 and mx!1 vectors, 
respectively. Each column of A is denoted A, and each scalar element of vector x is 


denoted x,, for i=1,...,n. Then, the LS problem has the equivalent form 
é€=min Ax, + A,X, +...+A,x, —b) (85) 


Suppose that x is an approximation to the solution x° after k iterations and its 
elements are x), for i=1,...,n. Considering an arbitrary element x,, all elements 


X,,.--X,, are updated in k+1 iterations while the rest x,,...x, are updated in k 


n 


iterations. Now, (85) can be written as 


n 
(k+l) (k+l) (k) 
MO 4 AO 4 DY AXP —b (86) 
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for the (k+1)" iteration of x“*”, which gives the local solution for x“*”. The 


argument of the term on the right hand side of (86) can be written as 


n 
(k+1) (k+1) (k) = 
Apxy + Axe + DY Ax} —b= 


Mz 


j=l 3 j=i+l ; (87) 
A — x?) + DIAS + YAS? —b 
j=l jzi 
By substituting s“? as the step or correction 
5) - xt — x) (88) 
(k,i-1) ; 
and r as the residual 
bid Sy aay Xe (k) 
i-l) _ + 
re = Aa + VARY —b (89) 
j=l jzi 
into (86) yields 
: (bk), (ki-1) 
é= min|A,s\ +r : (90) 





Thus, the global problem of finding the solution for x,,...,x, in (85) is equivalent 


n 


to solving the subproblems of (90), which can be assigned to the corresponding sensor 


nodes. Indeed, A, is known locally to the sensor node i, s\“ 


L 


’ is the locally computed 
correction, and r““"” is the residual after the i-1" node has completed its update x“, 
p p i-l 


The residual ee” is not locally available in the node i; however, it can be transmitted 


h 


by the i-1” sensor node. After the i‘" node calculates the new solution x‘"’, it sends 


the updated residual ro) to the other nodes. The residual can be shown to satisfy the 
p 


recursive equation [30] 


pi = ED So) (91) 


L 





whereas the new approximate solution is 


xP 1 45. (92) 
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Therefore, the i‘" node is assigned a column A, of the matrix A and assumes an 
initial solution x‘. Assuming that r, the initial estimation for the residual is available, 
the i‘" node solves (90) for s‘ and then updates its solution x‘ by (92). Following 
this step, the updated residual r”” is sent to the i+1” node in order to update its solution 


(1) 


x,,,. The procedure continues until a convergence criterion is satisfied. This series of 


approximations converge to the solution x’, and the norm of the residual decreases 


continuously [30], [31]. This procedure is summarized in Figure 25. 


a) Divide A intoitscolumns 4A,, for i=1,...,n, and assign one to each node 7. 
b) Initialize a? , for i=1,...,n, and ae 
c) For k=1 until convergence 

For i=1,...,n 


(k,i-1) |? 


Solve for s\ : minjA,sy? +r 
Sj 


Update the residual : r? =r“? + As. 


Send the updated residual to all nodes. 
Update the solution x@*? =x +5, 


i 


Check for convergence. 





Figure 25. Procedure for the distributed iterative solution of the LS problem (After 
Ref. [30]). 


This algorithm is highly distributed, and it can be implemented in a sensor 
network in order to spread the computational effort equally among all participating sensor 


nodes. 
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From the notation of (84), the matrix Ais the steering matrix D(0@)", each 
column A, is the i” column of D(@)", and x, is the weight coefficient w,. From (61), 


each column of D(@)" depends only on the position of the node and the set of 


approximation points. This column is available locally to the sensor node. Starting from 


an initial estimation for the residual r“” , the first node can solve (90) for s\ and update 


its own weight w<” by (92). Then, it can update the residual r’” and send it to the next 


node. However, since the residual is an mx1 vector, each transmission of the residual 
needs the transmission of m elements, which is a considerable amount of communication 
load. Another approach for the transmission of the residual seems more efficient. Assume 
that all sensors’ positions are initially broadcast to all sensor nodes so that each node 
knows the position of each other node; in this way, they can construct any of the columns 


th 


of the steering matrix. Then the 7” sensor node does not need to send the entire residual 


(k 


r“” but only the scalar correction s, and all other nodes can reconstruct the residual 


by repeating the (91). For example, the second node is not required to send the residual 
r“”) but just the scalar correction s, and then the third node will reconstruct the 


residual by computing (91): 
ae 2s (93) 


Therefore, in each iteration step, only one scalar element has to be broadcast, and 
the rest of the computations are performed locally. This distributed determination of the 


weights in a sensor network is summarized in Figure 26. 
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a) All nodes broadcast their position; each node can construct all columns 


D,, for i=1,...,.n, of D(O)". 


(0) 


b) Initialize w, i=1,...,n and r. 
c) For k =1 until convergence 
For i =1,...,n 
Receive s\ from i-1 node. 


7 k,i-l k,i-2 
Reconstruct the residual: r“"? =r" + D, 8. 


‘ 2 
. k,i-l 
Solve for s\” : min Degg el 
5; 


Update the residual: r“? =r" + Djs. 





Send the only the correction s\” to all nodes. 
Update the solution w*? = w +5. 


Check for convergence. 





Figure 26. Procedure for the proposed distributed iterative solution of the LS problem 
in a WSN environment. 


2. Computational and Communication Costs 


The performance of the distributed iterative procedure for the computation of the 
weights in a random sensor array is shown in Figures 27 and 28. In Figure 27, the 
residual norm is plotted for each local iteration for a 10x1 array of sensors. It is obvious 
that the residual norm converges to the actual residual after about thirty local iterations or 
about three complete iterations; one complete iteration is finished when all ten nodes 


perform a local iteration. Similar results are plotted in Figure 28 where the norm of the 
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(k) 


error e(k) between the approximate w and the actual solution w’ calculated after 


each complete iteration k is 





























(94) 
(e_ «-  --, —e Residual norm after approximation] _ 
-*- Residual norm for actual solution 
E ! : 
Z, i} | 
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Figure 27. Convergence of the residual norm to the actual residual, indicating that the 


algorithm converges to the real solution. After 3 complete iterations (30 local) the 
residual has converged. 
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Figure 28. Convergence of the norm of the error e(k) between the approximate w” 


and actual solution w’ calculated after each complete iteration k , indicating that 
the algorithm converges to the real solution. 


The computational cost can be measured in terms of the number of instructions 
needed for the implementation. Assuming that each node solves its local LS problem with 


a QR factorization obtained by Householder transformations, the cost per sensor NV, is 
given by (71), which for n = l yields 
N,, =2(m—-1/3)+m+1. (95) 


Note that the last two terms stand for the backsubstitution, so they are performed in each 


iteration, thus (95) becomes 


N,, = 2(m-1/3)+k(m+1). (96) 


Furthermore, each node has to reconstruct the mx1residual vector, so an additional 
number of 2m instructions are needed. Finally, the number of instructions per sensor 


results in 
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N 5 =2(m-1/3)+k(3m +1). (97) 
Therefore, the total processing cost for all the nodes N, is equal to 


N, =nxN,, =n[2(m-1/3)+k3m +1] (98) 


while the total processing power is 
P,=N,xP.=n[2(m-1/3)+k(3m+1)]x P. (99) 


Using the above equation, the processing cost in number of instructions for the 
implementation of the algorithm is plotted as a function of the number of sensors n, the 
number of approximation angles m, and the number of iterations k. In Figure 29, the 
number of instructions for the distributed approach become fewer than the centralized 
approach after a certain point, which is m=15; the simulation scenario is for a 10x1 
array of sensors and the cost has been computed for k = 5 iterations. Similarly, in Figure 
30, the processing cost for the centralized approach grows larger than the distributed 


approach as the number of the sensors increase (m= 20 and k =5). 
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Figure 29. Processing cost of the distributed algorithm as a function of the number of 
approximation angles for n=10 sensors and k =Siterations. Multiplying the 
number of instructions by P. gives the required processing power. 
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Figure 30. Processing cost of the distributed algorithm as a function of the number of 
sensors for m=20 approximation angles and k =5 iterations. Multiplying the 
number of instructions by P. gives the required processing power. 
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Figure 31. Processing cost of the distributed algorithm as a function of the number of 
iterations for n=10 sensors and m=20 approximation angles. Multiplying the 
number of instructions by P. gives the required processing power. 


Finally, Figure 31 shows the increase of the processing cost with the number of 
iterations, which indicates that the cost for the distributed approach grows higher than the 
cost for the centralized one after a specific number of iterations (k =5 in this case); a 


10x1 sensor array is deployed and twenty approximation points are used. 

The communication cost for the distributed approach can be derived as follows. 
Initially, each node broadcasts its own position, so a total of n data elements are 
transmitted. Then, after a local iteration step is finished, the scalar correction — has to 


be sent, so for k complete iterations, kn elements are transmitted. Therefore, the number 


of transmitted elements is 
C=(k+In (100) 


and it does not depend on the number of approximation angles m. Assuming that each 


element is represented by b bits, the total transmission power P. is given by 


Cc 
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P=CxP,xB=(k+l)nxP,xb. (101) 


In Figure 32, the communication cost is plotted as a function of the number of 
sensors. The increased communication load for the distributed approach is a reasonable 
tradeoff considering that this distributed approach relieves the central processing node 
from the entire computational load. Similarly, Figure 33 shows the effect of the number 
of iterations on the communication cost. Obviously, the speed of convergence positively 


affects the reduction of the transmitted elements. 
The total power for the implementation of the algorithm is the sum of the 
processing power given by (99) and the transmission power given by (101): 


P =n[2(m-1/3)+k(3m+1)]x P +(k+)nx P, xb (102) 


Considering the ratio7,, and substituting P, from (21), the normalized power P, 


as a number of required units of power P. is given by 


U 


P. = == n[2(m—1/3) + kGm-+ 1) P +(k+1)nxn, xb. (103) 


i 
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Figure 32. Communication Cost of the distributed algorithm as a function of the 
number of sensors for m=20 approximation angles and k =S iterations. 
Multiplying the number of data elements with P,xb gives the required 


transmission power. 
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Figure 33. Communication cost of the distributed algorithm as a function of the 
number of iterations for n=10 sensors and m=20 approximation angles. 
Multiplying the number of data elements with P,xb gives the required 
transmission power. 
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In Figures 34-36, the total normalized power needed for the implementation of the 
algorithm is plotted as a function of the number of approximation angles, the number of 
sensors and the number of iterations, respectively. In all cases, the power consumption of 
the distributed algorithm is larger than that of the centralized approach. However, the 
important characteristic, as in the previous algorithm, is that the total amount of power is 
shared among the nodes; consequently, the cluster head is relieved of the computational 
burden. For instance, in Figure 35, where m = 20 and k =5, the power for the distributed 
approach is almost three times larger than the centralized approach, but this is shared by 
n=20nodes. Thus, the average power consumption of an arbitrary node in the 
distributed network is six times less compared to the cluster head’s consumption in the 
centralized architecture. In a centralized approach, there is a considerably higher 
possibility that the cluster head will exhaust its limited battery and then the beamforming 
problem will have to be computed from scratch. In summary, the increased power 
consumption is reasonable enough to consider the distributed approach a viable solution 


if robustness is required in the network. 
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Figure 34. | Normalized power P, (number of P.) for both distributed and centralized 


approaches as a function of the number of approximation angles for 
n=10sensors and k = S iterations. 7,, is assumed to be 200 and b = 32 bits. 
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Figure 35. | Normalized power P, (number of P.) for both distributed and centralized 


approaches as a function of the number of sensors for m= 20 approximation 
points and k = 5 iterations. 77,, is assumed to be 200 and b = 32 bits. 
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Figure 36. | Normalized power P, (number of P.) for both distributed and centralized 


approaches as a function of the number iterations for n=10sensors and 
m = 20 approximation points. 77,, is assumed to be 200 and b = 32 bits. 


77 


C. SUMMARY 


In this chapter, two distributed approaches were proposed for the solution of the 
LS problem of the array weight computation. The presented algorithms were evaluated in 
terms of the processing cost, the communication cost and the total power consumption 


and then compared to the centralized implementation. 


In the first approach (distributed QR decomposition), the processing cost is the 
same as in the centralized approach, but there is a tradeoff of increased communication 
effort. However, only the first phase of the algorithm is power demanding, and potential 
modifications to the desired response do not require much power for the recalculation of 


the weights. 


In the second approach (iterative solution), there is rapid convergence to the 
actual solution, which yields a reduction of the processing cost if the number of iterations 
is kept low. The simulation results show that only 3 or 4 iterations are needed for the 
convergence of the algorithm, which results in considerably lower processing power. 
However, the communication cost is still higher when compared to the centralized 


approach. 


To sum it up, these two approaches require significantly lower average power per 


sensor node by efficiently sharing the power consumption among the nodes. 
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Vv. CONCLUSIONS 


The operational scenario adopted in this work assumes that a number of sensor 
nodes are randomly deployed in an area of interest in order to collect information about 
various kinds of objects. The acquired data has to be collected by an oveflying UAV; 
however, single sensors do not have sufficient power capabilities in order to establish 
communication with the UAV. Therefore, they are organized into clusters and cooperate 


in order to function as an array of sensor nodes. 


The effect of position errors on the performance of the random sensor array was 
analyzed, and the need for beamforming techniques which effectively mitigate that effect 
was discussed. Since reliability and robustness in a sensor network environment are 
crucial, two distributed algorithms for beamforming that efficiently manage to share the 
processing load among the sensor nodes, compared to centralized approaches which 
assign the entire effort to a single node, were presented. A simulation model was created 
and implemented in the MATLAB environment to evaluate the performance of the 


proposed algorithms. 


A. SIGNIFICANT RESULTS 


The simulations showed that the sidelobes in the array response increase as a 
function of the “randomness” of the sensor array. Thus, as the mean deviation from the 
uniform array was increased, the average sidelobe magnitudes also increased. These 
results validate the theoretical results of random arrays found in the literature. Another 
important point is that for the LS beamformer, a subset of approximation points can yield 
almost the same solution as larger sets. This offers significant reduction to the required 


processing and transmission power, which are crucial in sensor networks. 


Based on the performance analysis, the two proposed distributed algorithms can 
effectively share the processing load among the nodes. The first, a distributed 
implementation of the QR decomposition, has the same processing cost as the centralized 


one. The second approach, based on an iterative method of computing the weight vector 
19 


in the LS sense, converges quickly to the actual solution and achieves reduction of the 


total processing cost compared to that of the centralized one. 


For both algorithms, the tradeoff is the increased transmission power, causing an 
overall increase in the total power consumption in the network. This total power, 
however, is shared among the sensor nodes; therefore, the average power needed by a 
sensor node in the distributed implementation is lower than the power needed by the 
cluster head in the centralized approach. Consequently, the network’s susceptibility to 


failures due to excessive power consumption is greatly reduced. 


B. FUTURE WORK 


Throughout this work, several assumptions were made, such as the nodes can 
compute their positions without errors, and the communication between them is not 
affected by noise. A future effort may examine the effect of these errors on the 


computation of the weight vector and consequently on the array performance. 


In this work, the set of approximation points were selected based on uniform 
sampling, but there are other choices, such as using a non-uniform grid, which may 
require fewer points with similar performance. Initial results showed that certain 
approximation points, which have the physical meaning of direction angles, may be more 
important for the approximation of the desired response than others. These issues may be 


further investigated in a future study. 


The topology of WSN changes dynamically due to frequent additions and 
withdrawals of sensor nodes; some of them may switch in or off sleep mode and some 
other may fail because of the harsh environmental conditions or because of exhausted 
battery. For these scenarios, the processing and communication cost for the update of the 
weight vectors can be investigated. Also for the distributed algorithms, there is no need to 
solve the beamforming problem from scratch; if the array topology changes, it is 
important to analyze the effects on the costs and consequently the power needed for 


modifying the weight vector after a node has added or withdrawn from the array. 
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In this work, the emphasis was given to the data independent beamforming 
techniques, such as the LS approximation of the desired response. There are many proven 
data dependent techniques for which the weight vector can be determined adaptively 


[11], [12] and a distributed implementation of these methods could be examined. 


An array of M elements can be used to form a beampattern with exactly M -1 
narrowband nulls. These nulls can be placed towards the directions of incoming 
interferences in order to suppress them. Since there are straightforward polynomial based 
techniques for placing those nulls where desired in the array space, it would be 


interesting to investigate them in a future work. 


Finally, the communication cost was defined as a function of the data elements 
that need to be transmitted for the implementation of the algorithms. However, there is 
also a networking cost which consists of parameters such as the packet overhead and the 
retransmissions due to collisions and errors, which all add to the power consumption. 
Since the implementation of a distributed beamforming algorithm may be prohibitive by 
a high networking cost, it will be interesting to investigate the effect of the networking 


cost to the overall power consumption. 
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APPENDIX. 


MATLAB SOURCE CODE 


This appendix lists all MATLAB programs used in this work 






































e Array2D.m: 
SSS Filename: Array2D.m 
SS%S Author: Nikolaos Papalexidis 
53% Hellenic Air Force 
53% Date: June 2007 
SSS Description: This file generates the array beampattern for an 
SSS array with randomly positioned elements 
65% The weights are computed using the LS approach of 
SSS the desired response 
clear a 
close a 
ele 
SESESESSCEESESS PARAMETERS SESEESEESSEESSESES 
global c f 1 b Im Nx Ny 
c=3e8; 
f=2e9; 
l=c/f; 
b=2*pi/1; 
GdBavg=zeros (181,181); 
GdBerravg=zeros (181,181); 
GdBerrlinavg=zeros (181,181); 
GdBrefavg=zeros (181,181); 
GdBlsavg=zeros (181,181); 
GdBLSavg=zeros (181,181); 
GdBiteravg=zeros (181,181,16); 
GdBLSlavg=zeros (181,181); 
GdBLS2avg=zeros (181,181); 
GdBLS3avg=zeros (181,181); 
GdBLS4avg=zeros (181,181); 
GdBLS5Savg=zeros (181,181); 
wwls=[]; 
wwnu=[]; 
wwli=[]; 
WwW2=[]; 
WW3=[]; 
ww4=[]; 
WW5=[]; 
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6556655 SSS Uniform array (reference) 6656556555 


dx=1/2; % ideal distance lamda/2 in x-direction 
xn=(0:Nx-1) *dx; 

xn=repmat (xn',1,Ny); 

xn=reshape (xn,Nx*Ny,1); 


dy=1/2; % ideal distance lamda/2 in y-direction 
yn=(0:Ny-1) *dy; 

yn=repmat (yn, Nx,1); 

yn=reshape (yn, Nx*Ny,1); 





$ SSS SSSSS End of uniform array $SSESSESSSES 


for NI=1:NumIter; 
3%%%%%%S% DISTANCE DEVIATION %%%%3%%%3% 


if XG==1; 
devx=xe*l* (rand (Nx*Ny,1)-0.5); & Random deviation from xe% 
-lamda/2 to xe% lamda/2 
x=xntdevx; % Real positions in x-direction 





devy=ye*l* (rand (Nx*Ny,1)-0.5); % Random deviation from - 
lamda/2 to lamda/2 
y=yntdevy; % Real positions in x-direction 





elseif XG==2;% 2nd option - Firstly construct x,y and then 
assume linear 
[x, y]=rand_inter_dist (Nx,Ny); 
end 





SSESSESS Estimated position with defined error %%%%%%% 
xerr=x.%* (l+xest* (2*rand(Nx*Ny,1)-1)); % error position in x- 
direction with respect to actual 

yerr=y.* (l+yest* (2*rand(Nx*Ny,1)-1)); % error position in y- 


direction with respect to actual 


$SSSSSS% END OF DISTANCE DEVIATION $SSSSSS6SSESS 
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Amplitudes 


ooo 
$ 3% 


& Amplitudes 


ones (Nx*Ny, 1) ; 


Im= 


correct weights 


. 2 
1 oO 


weights2 (x, y,theta0,phi0) 


wm 


deviation 


wrong weights, 


% 


, 


weights2 (xerr, yerr,theta0,phi0) 


from actual position 


wmerrlin 


wmerr 


xeS 


assume 


7 % wrong weights, 


weights2 (xn, yn, theta0, phi0) 


perfect linear 


(uniform array) 


Reference weights 
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weights2 (xn, yn, theta0,phi0) ; 
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Gain from correct weights 
Gain from wrong weights, 


% 


gain2D (wm,x,y); 


G= 


deviation from 


% 


gain2D (wmerr,x,y); 
actual position 


Gerrlin 


Gerr 


assumed 


Gain from wrong weights, 


% 


i 


gain2D (wmerrlin, x,y) 


perfect linear 


gain2D (wref,xn,yn); 


Gref 


Gain from correct weights (dB) 


% 


Gain from wrong weigths, 


% 


(dB) 


10*1log10 (Gerrlin/max (max (G) )); 


Gain from wrong weights, 


% 





(dB) 


(uniform 


Gain for reference 


% 





10*1lo0g10 (Gerr/max (max(G))); 


10*1log10 (Gref/max (max (G))); 
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BEZEL EEE OEE OE COO OOO OE OO OO COO OEE EE EES LESSEE SELES EE ESES 
SB 

%%%%%%%% LS ESTIMATION OF THE WEIGHTS, GIVEN DESIRED RESPONSE 
SEZESEELBSS 

SEEEEEE EEO EEE EE EE EEE EEE EEE EE EEE ELSES EELS EEE ELLE SELES 
SE 


theta=—-90:90; 
th=theta*pi/180; 


dn=exp (j*b* (xn*sin (th) *cos (phi0) tyn*sin(th) *sin(phiO))); % 
steering vector for ULA 


Fdes=wref'*dn; 


d=exp (j*b* (x*sin(th) *cos (phi0)+y*sin(th)*sin(phi0))); % 
steering vector 


ww=Fdes/d; 
ww=ww'; & or ww=inv(d*d') *d*Fdes'; 


wwls=[wwls ww]; 


Gls=gain2D (ww,x,y); 
GdB1ls=10*1lo0g10 (Gls/max (max (G) )); 


thO=P (8); 
$%%%6%S6%S LS with fewer approximation points %%%% 


6% Uniform spacing %% 


5% (1) %% 

dt=2; 
rla=th0:-dt:-90; 
rla=flipdim(rla,2); 
rlb=th0+dt:dt:90; 
rl=[rla rlb]; 


tl=rl*pi/180; 
DN1=exp (j*b* (xn*sin(t1) *cos (phi0) tyn*sin(t1l)*sin(phiO))); 


FDES1l=wref'*DN1; 





Dl=exp (j*b* (x*sin(t1) *cos (phi0) +y*sin(t1) *sin(phi0O))); 





wwl=FDES1/D1; 
wwl=wwl'; 


WW1=[WW1 wwlj; 
87 


GLS1l=gain2D (wwl,x,y); 
GdBLS1=10*1log10 (GLS1/max (max (G))); 


5% (2) %% 
dt=4; 


r2a=th0:-dt:-90; 

r2a=flipdim(r2a,2); 

r2b=th0tdt:dt:90; 

r2=[r2a r2b]; 

t2=r2*pi/180; 

DN2=exp (j*b* (xn*sin(t2) *cos (phi0) tyn*sin(t2)*sin(phiO))); 


FDES2=wref'*DN2; 





D2=exp (j*b* (x*sin(t2) *cos (phi0) +y*sin(t2) *sin(phiO))); 


ww2=FDES2/D2; 
ww2=ww2'; 





WW2=[WW2 ww2]; 


GLS2=gain2D (ww2,x,y); 
GdBLS2=10*1o0g10 (GLS2/max (max (G))); 


5% (3) %% 

dt=6; 
r3a=th0:-dt:-90; 
r3a=flipdim(r3a,2); 
r3b=th0+dt:dt:90; 
r3=[r3a r3b]; 


t3=r3*pi/180; 


DN3=exp (j*b* (xn*sin(t3) *cos (phi0) tyn*sin(t3)*sin(phiO))); 





FDES3=wref'*DN3; 


D3=exp (j*b* (x*sin(t3) *cos (phi0) ty*sin(t3) *sin(phiO))); 


ww3=FDES3/D3; 
ww3=ww3'; 





WW3=[WW3 ww3]; 


GLS3=gain2D (ww3,x,y); 
GdBLS3=10*10g10 (GLS3/max (max (G))); 
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5% (4) %% 

dt=8; 
r4a=th0:-dt:-90; 
r4a=flipdim(r4a,2); 
r4b=th0+dt:dt:90; 
r4=[r4a r4b]; 


t4=r4*pi/180; 
DN4=exp (j*b* (xn*sin(t4) *cos (phi0) tyn*sin(t4) *sin(phiO))); 


FDES4=wref'*DN4; 





D4=exp (j*b* (x*sin(t4) *cos (phi0) +ty*sin(t4) *sin(phiO))); 


ww4=FDES4/D4; 
ww4=ww4'; 





WW4=[WW4 ww4]; 


GLS4=gain2D (ww4,x,y); 
GdBLS4=10*1lo0g10 (GLS4/max (max (G))); 


%% (5) %% 
dt=10; 


r5a=th0:-dt:-90; 

roa=flipdim(r5a,2); 

rob=th0tdt:dt:90; 

r5=[r5a r5b]; 

t5=r5*pi/18s0; 

DN5=exp (j*b* (xn*sin(t5) *cos (phi0) tyn*sin(t5)*sin(phiO))); 


FDESS=wref'*DN5; 





D5=exp (j*b* (x*sin(t5) *cos (phiO0) +y*sin(t5) *sin(phiO))); 


ww5=FDES5/D5; 
ww5=ww5'; 





WW5=[WW5 ww5]; 


GLS5=gain2D (ww5,x,y); 
GdBLS5=10*1o0g10 (GLS5/max (max (G))); 


%%$% Non uniform %%%% 
thetal=th0-12:4:th0+12; 
theta2=—-90:15:90; 
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ft=(theta2<th0-12) | (theta2>th0+12); 
thetaf=[thetal theta2(ft)]; 
THETA=sort (thetaf) ; 





TH=THETA*pi/180; 





DN=exp (j*b* (xn*sin (TH) *cos (phi0) tyn*sin (TH) *sin(phiO))); 





F 


O 


ES=wref'*DN; 





D=exp (j*b* (x*sin (TH) *cos (phi0) +y*sin(TH) *sin(phi0))); 


WW=FDES/D; 
WW=ww'; 





wwnu=[wwnu WW]; 


GLS=gain2D (WW,x,y); 
GdBLS=10*10g10 (GLS/max (max (G) ) ); 


229 0 Oo o 


YOAWAAIAAAAAAAAAAAIIFAAI 7999979973799 9909999 9-999 9999-999 9-99- 9999-99999 9 9. 








plot (xn,yn,'o') 

hold on; 

plot (x,y, 'rx') 

hold on; 

plot (xerr,yerr, 'mp') 
grid on 


axis equal 
title('Fig.1 Sensor array','Fontsize',12); 
legend('Perfect linear', 'Actual Position', 'Wrongly estimated'); 


figure (20); 




















plot (theta, GdB(:,phi_ang+90+1), 'Linewidth',2); 

hold on; 

plot (theta, GdBerr(:,phi_ang+90+1),'r-.', 'Linewidth',2); 
hold on; 

plot (theta, GdBerrlin(:,phi_ang+90+1),'m:', 'Linewidth',2); 
hold on; 

plot (theta, GdBref (:,phi_ang+90+1),'g-', 'Linewidth',2); 
hold on; 

plot (theta, GdBls(:,phi_ang+90+1),'c-', 'Linewidth',2); 
axis([-85 85 -50 5]); 

grid on 

title('Fig.2 : Beampattern for N linear array elements and given 





\phi', 'Fontsize',12); 
90 


xlabel('\theta (degrees)','Fontsize',12); 
ylabel('Power Gain (dB),' ,'Fontsize',12); 
( 














legend('Correct', 'Wrongly estimated', 'Assumed perfect linear', 'Ideal 
linear','LS weights'); 


GdBavg=GdBavg+GdB; 
GdBerravg=GdBerravg+GdBerr; 
GdBerrlinavg=GdBerrlinavg+GdBerrlin; 
GdBrefavg=GdBrefavg+GdBref; 
GdBlsavg=GdBlsavg+GdBls; 
GdBLSavg=GdBLSavg+GdBLS; 











GdBLSlavg=GdBLSlavg+GdBLS1; 
GdBLS2avg=GdBLS2avg+GdBLS2; 
GdBLS3avg=GdBLS3avg+GdBLS3; 
GdBLS4avg=GdBLS4avg+GdBLS4; 
GdBLSSavg=GdBLS5avg+GdBLS5; 














end $6%%5%S End of loop for iterated computations (average Gain) 


o2.¢. 2 '9.'C..o. 


GdBavg=GdBavg/NumIter; 
GdBerravg=GdBerravg/NumIter; 
GdBerrlinavg=GdBerrlinavg/NumIter; 
GdBrefavg=GdBrefavg/NumIter; 
GdBlsavg=GdBlsavg/NumIter; 
GdBLSavg=GdBLSavg/NumIter; 


GdBLSlavg=GdBLSlavg/NumIter; 
GdBLS2avg=GdBLS2avg/NumIter; 
GdBLS3avg=GdBLS3avg/NumIter; 
GdBLS4avg=GdBLS4avg/NumIter; 
GdBLSS5avg=GdBLSS5avg/NumIter; 











figure (30); 

theta=—90:90; 

lot (theta, GdBavg(:,phi_ang+90+1), 'Linewidth',2); 
old on; 
lot (theta, GdBerravg(:,phi_ang+90+1),'r-.', 'Linewidth',2); 
old on; 
lot (theta, GdBerrlinavg(:,phi_ang+90+1),'m:','Linewidth',2); 
old on; 
lot (theta, GdBrefavg(:,phi_ang+90+1),'g-', 'Linewidth',2); 
old. en; 
plot (theta, GdBlsavg(:,phi_ang+90+1),'c-', 'Linewidth',2); 
grid on 
legend('Correct', 'Wrongly estimated', 'Assumed perfect linear', 'Ideal 
linear','LS weights'); 
axis([-85 85 -50 5]); 
title('Fig.3 : Average Beampattern for N linear array elements and 











BO: Jee "O: ta 8G: ar" 























given \phi', 'Fontsize',12); 
xlabel('\theta (degrees)','Fontsize',12); 
ylabel('Power Gain (dB),' ,'Fontsize',12); 
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figure (40); 
































theta=—90:90; 

plot (theta, GdBavg(:,phi_ang+90+1), 'Linewidth',2); 

hold on; 

plot (theta, GdBerravg(:,phi_ang+90+1),'r-.', 'Linewidth',2); 
hold on; 

plot (theta, GdBrefavg(:,phi_ang+90+1),'g-', 'Linewidth',2); 
hold on; 

plot (theta, GdBlsavg(:,phi_ang+90+1),'c-', 'Linewidth',2); 
axis([-85 85 -50 5]); 

title('Fig.4 : Average Beampattern for N linear array 
elements', 'Fontsize',12); 

xlabel('\theta (degrees)','Fontsize',12); 

ylabel('Power Gain (dB),' ,'Fontsize',12); 

grid on 


legend('Correct', 'Wrongly estimated','Ideal linear','LS weights"); 


figure (50); 






































theta=—90:90; 

plot (theta, GdBavg(:,phi_ang+90+1), 'Linewidth',2); 

hold on; 

plot (theta, GdBrefavg(:,phi_ang+90+1),'g-', 'Linewidth',2); 
hold on; 

plot (theta, GdBlsavg(:,phi_ang+90+1),'c-', 'Linewidth',2); 
hold on; 

plot (theta, GdBLSavg(:,phi_ang+90+1),'k:', 'Linewidth',2); 
axis([-85 85 -50 5]); 

title('Fig.5 : Average Beampattern for N linear array 
elements', 'Fontsize',12); 

xlabel('\theta (degrees)','Fontsize',12); 

ylabel('Power Gain (dB),' ,'Fontsize',12); 

grid on 


legend('Correct', 'Ideal linear','LS weights','LS less constraints"); 
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e = =6Array2D.m: 


SSS Filename: InputParameters.m 

SSS Author: Nikolaos Papalexidis 

SSS Hellenic Air Force 

SSS Date: June 2007 

SSS Description: This function creates a GUI for defining the 


ole 
ol? 
ole 


characteristics of a random array 


function answer=inputparameters; 


prompt={'Give Nx number of array elements in x direction:',... 
"Give Ny number of array elements in y direction:',... 
"Generation of position: (1):Deviation from perfect linear, (2): 
From SCcracch", ss 
"Position error (with respect to perfect linear) in x direction 








"Position error (with respect to perfect linear) in y direction 








"Estimated position error (with respect to actual) in x 
direction (%$)", s«. 

"Estimated position error (with respect to actual) in y 
CPPECELOM. (29) pox 4 

"Elevation angle (theta) (degrees):',... 


"Azimuth angle (phi) (degrees):',... 
"Angle phi for beampattern:',... 

"Angle error in theta (+- degrees)',... 
"Angle error in phi (+-degrees)',... 
"Number of iterations', }; 





name='Parameters for antenna array'; 

numlines=1; 

défaultanswer={ Tho, %1", ti", oo, t207 Ot 2? . 2 B6*, Tag", 246, to, fot, 82 
S'h; 


answer=inputdlg (prompt,name, numlines, defaultanswer) ; 


for i=1l:length (answer) ; 
temp (i)=str2num(answer{i}); 
end 


answer=temp; 
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e Gain2D.m: 


SSS Filename: Gain2D.m 

SSS Author: Nikolaos Papalexidis 

SSS Hellenic Air Force 

SSS Date: June 2007 

SSS Description: This function calculates the beampattern gain of 


ol? 
ole 
ol? 


an array 


function Gain=gain2 (wm, xm, ym) ; 
global b 


for theta=—-90:90; 
for phi=-90:90; 
th=theta*pi/180; 
ph=phi*pi/180; 








F(90+thetat1, 90+phi+1)=sum(sum(conj (wm) .*exp (j*b* (xm*sin (th) *cos (ph) +ym 
*sin(th)*sin(ph))))); 

end 
end 


Gain=abs(F).%2; 
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e weights2.m : 


SSS Filename: weights2.m 

SSS Author: Nikolaos Papalexidis 

SSS Hellenic Air Force 

SS%S Date: June 2007 

SSS Description: This function calculates the weights for a 
uniform 

SSS array 


function w=weights2 (xm, ym, Theta, Phi); 


global Im b 


oe 


w=Im. *exp (j*b* (xm*sin (Theta) *cos (Phi) tym*sin (Theta) *sin(Phi))); 
weights 
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e rand_inter_dist.m: 


SSS Filename: rand_inter_dist.m 

SSS Author: Nikolaos Papalexidis 

SSS Hellenic Air Force 

SSS Date: June 2007 

SSS Description: This function creates a random array. the 
deviations 

SSS from the ideal array follow a uniform 
distribution 


function [xx, yy]=rand_inter_dist (Nx,Ny) ; 
global 1 


Ny=Ny+1; & First line will be ignored 
NxX=Nx+1; 








xx=zeros (Ny, Nx); 
yy=zeros (Ny, Nx) ; 


for i=1:Ny; 
for j=1:Nx; 
if j== 
xX (1,4) =0; 
else 
xx (i,j) =xx (i, j-1)+(rand*1/2+1/4) ; 
end 
if i==1; 
yy (i, 3)=0; 
else 
yy (1,3) =yy (i-1, 3) + (rand*1/2+1/4) ; 
end 
end 
end 
if Ny~=1 


XX1l=xx(2:Ny,2:Nx); 
yyl=yy (2:Ny,2:Nx) ; 
end 


xxl=reshape (xxl, (Nx-1) * (Ny-1),1); % ignore first line 
yyl=reshape (yyl, (Nx-1) * (Ny-1),1); 


XX=XxX1; 
yy=yyl; 
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e CostAnalysis.m : 
SSS Filename: CostAnalysis.m 
SSS Author: Nikolaos Papalexidis 
SSS Hellenic Air Force 
SSS Date: June 2007 
SSS Description: This file is used for the calculation of the 
SSS processing and communication costs and the 
SSS power consumption for the 
SSS centralized,the distributed QR decomposition 
SSS and the iterative approach. 
clear all 
close all 
CLG 
N=10; S Number of sensors 
M=10:30; % Number of angles 
6% Reference 
6% Computational Cost 
Comp_ref=2*N*2* (M-N/3) +M*N+N%2; 
Com_ref1=2*N; 
CESECECCECEEEEEEEEEEEEEEEEEEEEEEEEEEE EES SESE SESE SESE SESEEEECEEEEEES 
$$%S%SSS%SSSS6SS6SS = Distributed OR decomposition %%%%%%%%%%%S%S%SS%S 
CESECCCEECEEEEEEEEEEEEEEEEEEEEEEEEEEE SEES ES ESSE ESSE SESESEEEEEEEEEEES 
6% COMPUTATIONAL COST 
6% Number of Operations) 
6% Note that there is no distinguish between the type of operations 
6% like additions or multiplications 
% A matrix (M x N) 
% b vector (M x 1) 
% QR decomposition : Cqr=2*N*2* (M-N/3) 
% Update >: Cu=M*N 
% Back Substitution: Cbh=N%*2 


Cqr=2*N*2* (M-N/3); 
Cu=M*N; 

Cbh=N%*%2; 

Cub=CutCb; 


Comp1l=Cqr+Cub; 


figure (1); 


semilogy (M,Compl1, 'bo-',M,Cqr,'rp-.',M,Cub, 'md--', 'Linewidth',1.5); 
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title('Fig.1 : Computational Cost as a function of M, 

(N=10)', 'Fontsize',12); 

xlabel('Number of angles in beampattern approximation', 'Fontsize',12); 
ylabel('Number of instructions','Fontsize',12); 

legend('Total','QR decomposition', 'Update and Back Substitution'); 
grid on; 














M=20; 

N=5:20; 

Cqr=2*N.*2.* (M-N/3); 
Cu=M*N; 

Cbh=N.%2; 

Cub=CutCb; 


Comp2=Cqr+Cub; 


figure (2); 

semilogy (N,Comp2, 'bo-',N,Cqr,'rp-.',N,Cub, 'md--', 'Linewidth',1.5); 

title('Fig.1 : Computational Cost as a function of M, 

(N=10)', 'Fontsize',12); 

xlabel('Number of sensors','Fontsize',12); 

ylabel('Number of instructions','Fontsize',12); 

legend('Total','QR decomposition', 'Update and Back Substitution'); 
n 














, 


%% COMMUNICATION COST 
6% Defined as the number of data values we need to send (broadcasting) 
6% Not the number of packets 


% lst pass : From lst sensor to the last (QR decomposition) 


% Cqr= (M+4-N/2) * (N-1) 

& 2nd pass : Back substitution phase 
% Chs=2* (N-1) 

N=10; S Number of sensors 


M=10:30; % Number of angles 


Comm_gqr=(M+2-N/2) * (N-1); 
Comm_bs=2* (N-1); 


Com1l=Comm_qr+Comm_bs; 


figure (10); 

semilogy (M,Com1, 'bo-',M,Comm_qr, 'rp- 

.',M,ones (length (M),1)*Comm_bs, 'md:',M, ones (length (M),1) *Com_refl, 'gs-— 
2", 'banewideh",1.5) ; 

title ('Fig.3 :Communication Cost as a function of M , 

(N=10)', 'Fontsize',12); 

xlabel('Number of angles in beampattern approximation', 'Fontsize',12); 
ylabel('Number of elements transmitted', 'Fontsize',12); 
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legend('Total','QR decomposition ', 'Update and Back 
Substitution', 'Centralized approach"); 
grid. on; 


M=20; 
N=5:20; 
Com_ref2=2*N; 


Comm_qr=(M+2-N/2) .* (N-1); 
Comm_bs=2* (N-1); 


Com2=Comm_qr+Comm_bs; 


figure(11); 

semilogy (N,Com2, 'bo-',N,Comm_gr, 'rp-.',N,Comm_bs, 'md:',N,Com_ref2, 'gs- 
2’, Tainewidtkh? ,.1:..5):5 

title ('*Fig.4 :Communication Cost as a function of N , 

(M=20)', 'Fontsize',12); 

xlabel('Number of sensors','Fontsize',12); 

ylabel('Number of elements transmitted', 'Fontsize',12); 
legend('Total','QR decomposition ', 'Update and Back 

Substitution', 'Centralized approach"); 

grid on; 














%%% Power analysis 
Nt p=200; 

Pi=1; 

Ptb=Ntp*Pi; 

B=32; 


fo) 


% total Power 


fe) 


6 vs M 
P1=Comp1*Pit+tCom1*B*Ptb; 
Prefl=Comp1*Pi+Com_ref1*B*Ptb; 





N=10; S Number of sensors 
M=10:30; %& Number of angles 


figure (20); 
semilogy (M,P1,'bo-',M,Prefl,'rp:', 'Linewidth',1.5); 
title('Fig.20 :Power analysis as a function of M, 
(N=10)', 'Fontsize',12); 
xlabel('Number of angles in beampattern approximation', 'Fontsize',12); 
ylabel('Power (number of Pi)','Fontsize',12); 
legend('Distributed approach', 'Centralized approach"); 
n 














a. 


P2=Comp2*Pit+tCom2*B*Ptb; 
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Pref2= 


M=20; 
N=5:20 


figure 
semilo 
title ( 
(M=20) 
xlabel 
ylabel 
legend 
grid o 














2. 9.2. 2-2. 2. 


2. 9.2. 2-2. 


CP1=N* 


figure 
semilo 
title ( 
(N=10, 
xlabel 
ylabel 
legend 
grid o 














CP2=N* 


figure 
semilo 
title ( 
(M=20) 
xlabel 
ylabel 
legend 














Comp2*Pi+Com_ref2*B*Ptb; 


i 


(25); 

gy (N,P2, 'bo-',N,Pref2,'rp:', 'Linewidth',1.5); 
"Fig.25 iPower analysis as a Function of N , 

', 'Fontsizge*,12); 

('Number of sensors', 'Fontsize',12); 

("Power (number of Pi)', 'Fontsize',12); 
("Distributed approach', 'Centralized approach'); 
n 


1, 


9°099099099099099099099099099099090909909909909090909090909090909090990909099090090090 
SECC SSCS SSE SSE SCE S CESSES CESSES CESSES CES CESSES CESS 
ao2o0o00009000 0o09009009009009009009000900 
$%%SSS%S6%S Iterative Approach %6SSSSSSSSSSSSSCSCESS 
o0o00000909090909090909090909090909090909090909090909090909090909090909090909090909090909999 
SES CS SSE SSE SSE SCE SSE SCE SCE SSCS CESCE CCE SCE ECESCE CSS 


%S Number of sensors 
0; % Number of angles 


(2* (M-1/3) +K* (3*M+1)); 


(30); 

gy (M,CP1, 'bo-',M,Comp1, 'rp:', 'Linewidth',1.5); 

"Fig.30 :Computational Cost as a function of M , 

K=5)', 'Fontsize',12); 

("Number of angles in beampattern approximation', 'Fontsize',12); 
(‘Number of instructions', 'Fontsize',12); 

("Distributed approach', 'Centralized approach"); 

ny; 


(2* (M-1/3) +K* (3*M+1)); 


(35); 

gy (N,CP2, 'bo-',N,Comp2, 'rp:', 'Linewidth',1.5); 
'Fig.35 :Computational Cost as a Function of N , 
‘, "Fontsize", 12); 

("Number of sensors', 'Fontsize',12); 

('Number of instructions','Fontsize',12); 
("Distributed approach', 'Centralized approach"); 


100 


CP3=N* (2* (M-1/3) +K* (3*M+1) ); 
Cref=2*N*2* (M-N/3) +M*N+N%2; 


figure (40); 

semilogy (K,CP3, 'bo-',K,Cref*ones (length(K),1),'rp:','Linewidth',1.5); 
title('Fig.40 :Computational Cost as a function of iterations K , 
(N=10,M=20)', 'Fontsize',12); 

xlabel('Number of iterations','Fontsize',12); 

ylabel('Number of instructions','Fontsize',12); 

legend('Distributed approach', 'Centralized approach"); 

grid on; 














N=10; S Number of sensors 
M=10:30; % Number of angles 
K=5; 


ComParl=(K+1) *N; 
figure (50) 

M=20 
N=5: 
K=5; 


, 


20; 


ComPar2=(K+1) *N; 


semilogy (N,ComPar2, 'bo-',N,Com_ref2,'rp:', 'Linewidth',1.5); 
title('Fig.50 :Communication Cost as a function of N , 
(M=20,K=5)', 'Fontsize',12); 

xlabel('Number of sensors', 'Fontsize',12); 

ylabel('Number of elements transmitted', 'Fontsize',12); 
legend('Distributed approach', 'Centralized approach"); 

grid on; 














figure (60); 


ComPar3=(K+1) *N; 
Com_ref3=2*N; 
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semilogy (K,ComPar3, 'bo- 

',K,ones (length (K),1)*Com_ref3,'rp:', 'Linewidth',1.5); 
title('Fig.60 :Communication Cost as a function of K , 
(N=10,M=20)', 'Fontsize',12); 

xlabel('Number of iterations','Fontsize',12); 
ylabel('Number of elements transmitted', 'Fontsize',12); 
legend('Distributed approach', 'Centralized approach"); 
grid on; 














66% Power Analysis 
Nt p=200; 

Pi=1; 

Ptb=Ntp*Pi; 

B=32; 


N=10; S Number of sensors 
M=10:30; & Number of angles 
K=4; 


Parl=CP1*Pit+ComPar1*B*Ptb; 
Prefl=Comp1*Pi+Com_ref1*B*Ptb; 


figure (70); 

plot (M, Parl, 'bo-',M,Prefl1,'rp:', 'Linewidth',1.5); 

title('Fig.70 :Power analysis as a function of M, 

(N=10)', 'Fontsize',12); 

xlabel('Number of angles in beampattern approximation’, 'Fontsize',12); 
ylabel('Power (number of Pi)','Fontsize',12); 

legend('Distributed approach', 'Centralized approach"); 

grid on; 














Par2=CP2*Pit+ComPar2*B*Ptb; 
Pref2=Comp2*Pi+Com_ref2*B*Ptb; 


figure (80); 

semilogy (N,Par2, 'bo-',N,Pref2,'rp:', 'Linewidth',1.5); 
title('Fig.80 :Power analysis as a function of N , 
(M=20)', 'Fontsize',12); 

xlabel('Number of sensors', 'Fontsize',12); 
ylabel('Power (number of Pi)','Fontsize',12); 
legend('Distributed approach', 'Centralized approach"); 
grid on; 














M=20; 
N=10; 
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Par3=CP3*Pit+tComPar3*B*Ptb; 


Cqr=2*N.*2.%* (M-N/3) ; 
Cu=M*N; 
Ch=N.%2; 
Cub=CutCb; 


Comp3=Cqr+Cub; 


Pref3=Comp3*Pit+tCom_ref3*B*Ptb; 


figure (90); 














plot (K, Par3, 'bo-',K, ones (length (K),1)*Pref3,'rp:', 'Linewidth',1.5); 
title('Fig.90 :Power analysis as a function of K , 

(M=20)', 'Fontsize',12); 

xlabel('Number of iterations','Fontsize',12); 

ylabel('Power (number of Pi)','Fontsize',12); 

legend('Distributed approach', 'Centralized approach"); 

grid on; 


103 


AJP oP oP oP 
JP oP oP oN? 
AJP AP oP oP 


ol? 
ole 
ole 


ol? 
ol? 
ol? 


e Iterative.m : 


Filename: 
Author: 


Date: 


Description: 


Iterative.m 

Nikolaos Papalexidis 

Hellenic Air Force 

June 2007 

This file is used for the implementation of the 
distributed iterative solution of the LS problem 





XPOS=x; 
YPpos=y,; 


A=D5'; 
q=FDES5'; 





[M,N]=size (A); 


Clear 21 22 23 24 25 26 2? 2@ 2z mn 4 


Z1=[]; 
Z2=[]; 
Z3=[]; 
Z4=[]; 
Z5=[]; 
Z6=[]; 
Z7=[]; 








x 
~“e 


x 
~“e 


x 
~e 


x 
oo™Ne 


x 
~‘e 


~ ON 
yer rervrvr vr vr rH 
NeoNe x 


> 
oO 
ll 
D> 
~ 
AAIDOBWNHPR 


zzl=zeros (N,1 


zzZz1(1)=z1; 


zZz2=zeros(N,1 


ZZ2(2)=z2; 


zZz3=zeros(N,1 


223 (3)=23; 








zz4=zeros (N,1 
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zzZ4(4)=24; 











ZZ5=zeros(N,1); 
220 (5) =257 
zz6=zeros(N,1); 
ZzZ6(6)=z6; 
zz/=zeros(N,1); 
ZzzZz7/7(7)=z7; 
zz8=zeros(N,1); 
z2z8(8)=z8; 
s1=0; 

s2=0; 

s3=0; 

s4=0; 

s5=0; 

s6=0; 

s7=0; 

s8=0; 








ZZ=Z2214+222+223+224+225+226+227+228; 


r0=A*zz(:,1)-q; 

Niter=6; & Number of iterations 
r=zeros (M,Niter*N); % residual 
r(:,1)=r0; 

rhat=r0; 

z=A\q; 


for k=1:Niter; 
m=(k-1) *8+1; 


s1=-Al\rhat; 
rhat=rhattAl*sl1; 
r(:,mtl)=rhat; 
zl=zl+sl1; 
Z1=[(Z1 zl]; 





s2=-A2\rhat; 
rhat=rhattA2*s2; 
r(:,mt+2)=rhat; 
Z2=Z2+82; 
Z2=[Z2 22]; 





s3=-A3\rhat; 
rhat=rhatt+A3*s3; 
r(:,mt+3)=rhat; 
Z3=2z3+s3; 
Z23=[2Z3 23]; 





s4=-A4\rhat; 
rhat=rhattA4*s4; 
r(:,mt4)=rhat; 
zZ4=z4+s4; 
Z4=[Z4 24]; 
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s5=-A5\rhat; 
rhat=rhat+A5*s5; 
r(:,mt+5)=rhat; 
Z5=z5+s5; 
Z25=[Z25 z5];; 





s6=-A6\rhat; 
rhat=rhatt+A6*s6; 
r(:,mt+6)=rhat; 
zZ6=z6+sS6; 
Z6=[Z6 z6]; 





s7=-A7\rhat; 
rhat=rhattA7*s7; 
r(:,mt+7)=rhat; 
zZ7=z7+s7; 
Z7=[Z7 z7]; 





s8=-A8\rhat; 
rhat=rhatt+A8*s8; 
r(:,m+8)=rhat; 
zZ8=z8+s8; 
Z8=[Z8 z8]; 








zz(:,k)=[21;22;23;24;25;2z6;27; 28]; 
end 


for i=1:Niter*N+1; 
n(i)=norm(r(:,i)); 
end 


R=norm(A*z-q) ; 


figure (1); 
tl=1l:length (21); 
Nl=length (tl); 


subplot (3,3,1); 

plot (t1,ones(N1,1)*real(z(1)),'r-',t1l,real(Z1),'b-."'); 
title('Real of w_1', 'Fontsize',12); 

grid on; 


subplot (3,3,2); 

plot (t1,ones(N1,1)*real(z(2)),'r-',t1l,real(Z2),'b-.'); 
title('Real of w_2', 'Fontsize',12); 

grid on; 


subplot (3,3,3); 

plot (t1,ones(N1,1)*real(z(3)),'r-',t1l,real(Z3),'b-.'); 
title('Real of w_3', 'Fontsize',12); 

grid on; 
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subplot (3,3,4); 

plot (t1,ones(N1,1)*real(z(4)),'r-',tl,real(Z4),'b-.'); 
title('Real of w_4', 'Fontsize',12); 

grid on; 


subplot (3,3,5); 

plot (t1,ones(N1,1)*real(z(5)),'r-',tl,real(Z5),'b-.'); 
title('Real of w_5', 'Fontsize',12); 

grid-on; 


subplot (3,3,6); 

plot (t1,ones(N1,1)*real(z(6)),'r-',t1l,real(Z6),'b-.'); 
title('Real of w_6', 'Fontsize',12); 

grid on; 


subplot (3,3,7); 

plot (t1,ones(N1,1)*real(z(7)),'r-',t1l,real(Z7),'b-.'); 
title ('Real of w_7", "Fontsize',12); 

grid on; 


subplot (3,3,8); 

plot (t1,ones(N1,1)*real(z(8)),'r-',t1l,real(Z8),'b-.'); 
title('Real of w_8', 'Fontsize',12); 

grid on; 


subplot (3,3,9); 

plot (n); 

hold on; 

plot (ones (Niter*N+1,1)*R,'r-'); 
title('Residual Norm', 'Fontsize',12); 
grid. on; 





figure (2); 
tl=1l:length (21); 
Nl=length (tl); 


subplot (3,3,1); 

plot (t1,ones(N1,1)*imag(z(1)),'r-',tl,imag(Z1),'b-.'); 
title('Imaginary of w_1','Fontsize',12); 

grid on; 


subplot (3,3,2); 

plot (t1,ones(N1,1)*imag(z(2)),'r-',tl,imag(Z2),'b-.'); 
title('Imaginary of w_2','Fontsize',12); 

grid en; 


subplot (3,3,3); 

plot (t1,ones(N1,1)*imag(z(3)),'r-',tl,imag(Z3),'b-.'); 
title('Imaginary of w_3','Fontsize',12); 

grid on; 
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subplot (3,3,4); 


plot (t1,ones(N1,1)*imag(z(4)),'r-',tl,imag(Z4), 


title('Imaginary of w_4','Fontsize',12); 
grid on; 


subplot (3,3,5); 


plot (t1,ones(N1,1)*imag(z(5)),'r-',tl,imag(Z5), 


title('Imaginary of w_5', 'Fontsize',12); 
grid-on; 


subplot (3,3,6); 


plot (t1,ones(N1,1)*imag(z(6)),'r-',tl,imag(Z6), 


title('Imaginary of w_6','Fontsize',12); 
grid on; 


subplot (3,3,7); 


plot (t1,ones(N1,1)*imag(z(7)),'r-',tl,imag(Z7), 


title('Imaginary of w_7','Fontsize',12); 
grid on; 


subplot (3,3,8); 


plot (t1,ones(N1,1)*imag(z(8)),'r-',tl,imag(Z8), 


title('Imaginary of w_8', 'Fontsize',12); 
grid-on; 


subplot (3,3,9); 

plot (n); 

hold on; 

plot (ones (Niter*N+1,1)*R,'r-'); 
title('Residual Norm', 'Fontsize',12); 
grid on; 





figure (3); 
tl=1l:length (21); 
Nl=length (tl); 


subplot (3,3,1); 


plot (t1,ones(N1,1)*abs(z(1)),'r-',tl,abs(Z21), 


title('Magnitude of w_1','Fontsize',12); 
grid on; 


subplot (3,3,2); 


plot (t1,ones(N1,1)*abs(z(2)),'r-',tl,abs(Z2), 


title ('Magnitude of w_2','Fontsize',12); 
grid on; 


subplot (3,3,3); 


plot (t1,ones(N1,1)*abs(z(3)),'r-',t1l,abs(Z3), 


title ('Magnitude of w_3','Fontsize',12); 
grid on; 
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subplot (3,3,4); 

plot (t1,ones(N1,1)*abs(z(4)),'r-',tl,abs(Z4),'b-.'); 
title ('Magnitude of w_4','Fontsize',12); 

grid on; 


subplot (3,3,5); 

plot (t1,ones(N1,1)*abs(z(5)),'r-',tl,abs(Z5),'b-.'); 
title ('Magnitude of w_5', 'Fontsize',12); 

grid on; 


subplot (3,3,6); 

plot (t1,ones(N1,1)*abs(z(6)),'r-',t1l,abs(Z6),'b-.'); 
title ('Magnitude of w_6','Fontsize',12); 

grid on; 


subplot (3,3,7); 

plot (t1,ones(N1,1)*abs(z(7)),'r-',tl,abs(Z7),'b-.'); 
title ('Magnitude of w_7','Fontsize',12); 

grid on; 


subplot (3,3,8); 

plot (t1,ones(N1,1)*abs(z(8)),'r-',t1l,abs(Z8),'b-.'); 
title ('Magnitude of w_8', 'Fontsize',12); 

grid on; 


subplot (3,3,9); 

plot (n); 

hold on; 

plot (ones (Niter*N+1,1)*R,'r-'); 
title('Residual Norm', 'Fontsize',12); 
grid. on; 





figure (4); 

semilogy(n,'bo-'); 

hold on; 

semilogy (ones (Niter*N+1,1)*R,'rp-'); 
title('Residual Norm', 'Fontsize',12); 
xlabel('Number of local iterations','Fontsize',12); 
ylabel('Residual Norm', 'Fontsize',12); 

grid on; 




















figure (5); 


for i=1:Niter; 
ztemp=[21(i);22(i);23(1);24(i);25(1);26(1);27(1);28 (1) ]; 
e(1)=norm(z-ztemp) ; 

end 





figure (6) 
semilogy(e,'bo-'); 
title('Norm of the weight error','Fontsize',12); 
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x] 
y] 


abel 








abel 








("Number of complete iterations', 'Fontsize',12); 
("Norm of the error between approximate and actual 


weights', 'Fontsize',12); 
grid on; 
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